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Abstract. Unlike equilibrium statistical mechanics, with its well-established 
foundations, a similar widely- accepted framework for non-equilibrium statistical 
mechanics (NESM) remains elusive. Here, we review some of the many recent 
activities on NESM, focusing on some of the fundamental issues and general 
aspects. Using the language of stochastic Markov processes, we emphasize general 
properties of the evolution of configurational probabilities, as described by master 
equations. Of particular interest are systems in which the dynamics violate 
detailed balance, since such systems serve to model a wide variety of phenomena in 
nature. We next review two distinct approaches for investigating such problems. 
One approach focuses on models sufficiently simple to allow us to find exact, 
analytic, non-trivial results. We provide detailed mathematical analyses of a one- 
dimensional continuous-time lattice gas, the totally asymmetric exclusion process 
(TASEP). It is regarded as a paradigmatic model for NESM, much like the role 
the Ising model played for equilibrium statistical mechanics. It is also the starting 
point for the second approach, which attempts to include more realistic ingredients 
in order to be more applicable to systems in nature. Restricting ourselves to 
the area of biophysics and cellular biology, we review a number of models that 
are relevant for transport phenomena. Successes and limitations of these simple 
models are also highlighted. 



1. Introduction 

What can we expect of a system which consists of a large number of simple 
constituents and evolves according to relatively simple rules? To answer this question 
and bridge the micro-macro connection is the central goal of statistical mechanics. 
About a century ago, Boltzmann made considerable progress by proposing a bold 
hypothesis: When an isolated system eventually settles into a state of equilibrium, 
all its microstates are equally likely to occur over long periods of time. Known 
as the microcanonical ensemble, it provides the basis for computing averages of 
macroscopic observables: (a) by assuming (time independent) ensemble averages can 
replace time averages in such an equilibrium state, (b) by labeling each microstate C, 
(a configuration of the constituents which can be reached via the rules of evolution) 
as a member of this ensemble, and (c) by assigning the same weight to every member 
(P* (C) (X 1). This simple postulate forms the foundation for equilibrium statistical 
mechanics (EQSM), leads to other ensembles for systems in thermal equilibrium, and 
frames the treatment of thermodynamics in essentially all textbooks. The problem of 
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answering the question posed above shifts, for systems in equihbrium, to computing 
averages with Boltzmann weights. 

By contrast, there is no similar stepping stone for non-equihbrium statistical 
systems (NESM), especially ones far from equilibrium. Of course, given a set of rules 
of stochastic evolution, it is possible to write down equations which govern the time 
dependent weights, P{C,t). But that is just the starting point of NESM, as little is 
known, in general, about the solutions of such equations. Even if we have a solution, 
there is an added complication: the obvious inequivalence of time- and ensemble- 
averages a la Boltzmann. Instead, since our interest is the full dynamic behavior of 
such a statistical system, we must imagine (a) performing the same experiment many 
times, (b) collecting the data to form an ensemble of trajectories through configuration 
space, and (c) computing time dependent averages of macroscopic observables from 
this ensemble. The results can then be compared to averages obtained from P{C,t). 
Despite these daunting tasks, there are many studies [TJ [H with the goal of 
understanding such far-from-equilibrium phenomena. 

Here, we focus on another aspect of NESM, namely, systems which evolve 
according rules that violate detailed balance. In general, much less is known about 
their behavior, though they are used to model a much wider range of natural 
phenomena. Examples include the topic in this review - transport in biological 
systems, as well as epidemic spreading, pedestrian/vehicular traffic, stock markets, 
and social networks. A major difficulty with such systems is that, even if such a system 
is known (or assumed) to settle eventually in a time-independent state, the appropriate 
stationary weights are not generally known. In other words, there is no overarching 
principle, in the spirit of Boltzmann's fundamental hypothesis, which provides the 
weights for such non- equilibrium steady states (NESS). We should emphasize that, if 
the dynamics is Markovian, then these weights can be constructed formally from the 
rules of evolution [4]. However, this formal solution is typically far too intractable to 
be of practical use. As a result, such NESS distributions are explicitly known only 
for a handful of model systems. Indeed, developing a fundamental and comprehensive 
understanding of physics far from equilibrium is recognized to be one of the 'grand 
challenges' of our time, by both the US National Academy of Sciences [S] and the US 
Department of Energy [6]. Furthermore, these studies point out the importance of 
non-equilibrium systems and their impact far beyond physics, including areas such as 
computer science, biology, public health, civil infrastructure, sociology, and finance. 

One of the aims of this review is to provide a framework in which issues of NESM 
are well-posed, so that readers can appreciate why NESM is so challenging. Another 
goal is to show that initial steps in this long journey have been taken in the form 
of a few mathematically tractable models. A good example is the totally asymmetric 
simple exclusion process (TASEP). Like the Ising model, TASEP also consists of binary 
constituents, but evolves with even simpler rules. Unlike the scorn Ising's model 
faced in the 1920's, the TASEP already enjoys the status of a paradigmatic model. 
Fortunately, it is now recognized that seemingly simplistic models can play key roles 
in the understanding fundamental statistical mechanics and in formulating applied 
models of real physical systems. In this spirit, our final aim is to provide potential 
applications of the TASEP, and its many relatives, to a small class of problems in 
biology, namely, transport at molecular and cellular levels. 

This article is organized along the lines of these three goals. The phrase 'non- 
equilibrium statistical mechanics' has been used in many contexts, referring to very 
different issues, in a wide range of settings. The first part of the next section will help 
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65 readers discern the many facets of NESM. A more specific objective of section [3] is to 

66 review a proposal for characterizing ah stationary states by a pair of time independent 

67 distributions, {P* {C) ,K* {C ^ C')}, where K* is the probabihty current 'flowing' 

68 from C to C 7 . In this scheme, ordinary equilibrium stationary states (EQSS) 

69 appear as the restricted set {P*,0}, whereas states with K* ^ are identified as 

70 NESS. Making an analogy with electromagnetism, this distinction is comparable to 

71 that of electrostatics vs. magnetostatics, as the hallmark of the latter is the presence 

72 of steady and persistent currents. Others (e.g., [8]) have also called attention to 

73 the importance of such current loops for NESS and the key role they play in the 

74 understanding of fluctuations and dissipation. Two examples of NESS phenomena, 

75 which appear contrary to the conventional wisdom developed from EQSM, will be 

76 provided here. 

77 Section [3] will be devoted to some details on how to 'solve' the TASEP, for 

78 readers who are interested in getting involved in this type of study. In particular, 

79 we will present two complementary techniques, with one of them exploiting the 

80 relationship between two-dimensional systems in equilibrium and one-dimensional 

81 systems in NESM. While TASEP was introduced in 1970 [9] for studying interacting 

82 Markov processes, it gained wide attention two decades later in the statistical physics 

83 community |10l 111) I12j . In a twist of history, two years before its formal introduction, 

84 Gibbs and his collaborators introduced [T31 [13] a more complex version of TASEP to 

85 model mRNA translation in protein synthesis. As the need for modeling molecular 

86 transport in a biological setting provided the first incentives for considering such 

87 NESM systems, it is fitting that we devote the next part, section HI to potential 

88 applications for biological transport. In contrast to the late 60's, much more about 

89 molecular biology is known today, so that there is a large number of topics, even within 

90 this restricted class of biological systems. Though each of which deserves a full review, 

91 we will limit ourselves to a few paragraphs for each topic. The reader should regard 

92 our effort here as a bird's eye view 'tour guide', pointing to more detailed, in-depth 

93 coverages of specific avenues within this rich field. Finally, we should mention that 

94 TASEP naturally lends itself to applications in many other areas, e.g., traffic flow [TS] 

95 and surface growth [UllII], etc. All are very interesting, but clearly beyond the scope 

96 of this review, as each deserves a review of its own. In the last section, [SJ we conclude 

97 with a brief summary and outlook. 

98 2. General aspects of non-equilibrium statistical mechanics 

99 In any quantitative description of a system in nature, the first step is to specify the 

100 degrees of freedom to focus our attention, while ignoring all others. For example, in 

101 Galileo's study of the motion of balls dropped from a tower, the degrees of freedom 

102 associated with the planets is ignored. Similarly, the motion of the atoms within the 

103 balls plays no role. The importance of this simple observation is to recognize that 

104 all investigations are necessarily limited in scope and all quantitative predictions are 

105 approximations to some degree. Only by narrowing our focus to a limited window 

106 of length- and/or time-scales can we make reasonable progress towards quantitative 

107 understanding of natural phenomena. Thus, we must start by specifying a set 

108 of configurations, {C}, which accounts for all the relevant degrees of freedom of 

109 the system. For example, for the traditional kinetic theory of gases (in d spatial 

110 dimensions), C is a point in a 2dN dimensional phase space: {xi,pi \ , i = 1, N . For 

111 an Ising model with spins s = ±1, the set {C} is the 2^ vertices of an N dimensional 
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112 cube: {sj} , i = 1, N . In the first example, which should be suitable for describing 

113 Argon at standard temperature and pressure, the window of length scales certainly 

114 excludes Angstroms or less, since the electronic and hadronic degrees of freedom within 

115 an atom are ignored. Similarly, the window in the second example also excludes 

116 many details of solid state physics. Yet, the Ising model is remarkably successful at 

117 predicting the magnetic properties of several physical systems [18j [191 ED] ■ 

118 Now, as we shift our focus from microscopic to macroscopic lengths, both {C} 

119 and the description also change. Keeping detailed accounts of such changes is 

120 the key idea behind renormalization, the application of which led to the extremely 

121 successful prediction of non-analytic thermodynamic properties near second order 

122 phase transitions [21]. While certain aspects of these different levels of description 

123 change, other aspects - e.g., fundamental symmetries of the system - remain. One 

124 particular aspect of interest here is time reversal. Although physical laws at the 

125 atomic level respect this symmetrjfjj, 'effective Hamiltonians' and phenonienological 

126 descriptions at more macroscopic levels often do not. One hallmark of EQSM is that 

127 the dynamics, effective for whatever level of interest, retain this symmetry. Here, the 

128 concept and term 'detailed balance' is often used as well as 'time reversal.' By contrast, 

129 NESM provides a natural setting for us to appreciate the significance of this micro- 
no macro connection and the appearance of time-irreversible dynamics. We may start 

131 with a system with many degrees of freedom evolving with dynamics obeying detailed 

132 balance. Yet, when we focus on a subsystem with far fewer degrees of freedom, it is 

133 often reasonable to consider a dynamics that violates detailed balance. Examples of 

134 irreversible dynamics include simple friction in solid mechanics, resistance in electrical 

135 systems, and viscosity in fluid flows. 

136 Before presenting the framework we will use for discussing fundamental issues 

137 of NESM, let us briefly alert readers to the many settings where this term is used. 

138 Starting a statistical system in some initial conflguration, Co, and letting it evolve 

139 according to rules which respect detailed balance, it will eventually wind up in an 

140 EQSS (precise definitions and conditions to be given at the beginning of section [2] 

141 below). To be explicit, let us denote the probability to find the system in configuration 

142 C at time t by P{C,t) and start with P(C,0) = (5(C-Co). Then, P{C,t oo) 

143 will approach a stationary distribution, P*(C), which is recognized as a Boltzmann 

144 distribution in equilibrium physics. Before this 'eventuality', many scenarios are 

145 possible and all of them rightly deserve the term NESM. There are three important 

146 examples from the literature. Physical systems in which certain variables change so 

147 slowly that reaching P* may take many times the age of the universe. For time scales 

148 relevant to us, these systems are always 'far from equilibrium'. To study the statistics 

149 associated with fast variables, these slow ones might as well be considered frozen, 

150 leading to the concept of 'quenched disorder'. The techniques used to attack this class 

151 of problems are considerably more sophisticated than computing Boltzmann weights 

152 [3 [1], and are often termed NESM. At the other extreme, there is much interest in 

153 behavior of systems near equilibrium, for which perturbation theory around the EQSS 

154 is quite adequate. Linear response is the first step in such approaches [23l [24l [25l [26] . 

155 with a large body of well established results and many textbooks devoted to them. 

156 Between these extremes are system which evolve very slowly, yet tractably. Frequently, 

157 these studies come under the umbrella of NESM and are found with the term 'aging' 



I Strictly speaking, if we accept CPT as an exact symmetry, then time reversal is violated at the 
subatomic level, since CP violation has been observed. So far, there is no direct observation of T 
violation. For a recent overview, see, e.g., reference |22| . 
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158 in their titles 

159 Another frequently encountered situation is the presence of time-dependent 

160 rates. Such a problem corresponds to many experimental realizations in which 

161 control parameters, e.g., pressure or temperature, are varied according to some time- 

162 dependent protocol. In the context of theoretical investigations, such changes play 

163 central roles in the study of work theorems [27l[28l[Ml[3Ql[3ll[32l[33l|34l|35]. In 

164 general these problems are much less tractable and will not be considered here. 

165 By contrast, we will focus on systems which evolve according to dynamics that 

166 violate detailed balance. The simplest context for such a system is the coupling to two 

167 or more reservoirs (of the same resource, e.g., energy) in such a way that, when the 

168 system reaches a stationary state, there is a steady flux through it. A daily example 

169 is stove-top cooking, in which water in a pot gains energy from the burner and loses 

170 it to the room. At a steady simmer, the input balances the heat loss and our system 

171 reaches a non-equilibrium steady state (NESS). That these states differ significantly 

172 from ordinary EQSS's has been demonstrated in a variety of studies of simple model 

173 systems coupled to thermal reservoirs at two different temperatures. Another example, 

174 at the global scale, is life on earth, the existence of which depends on a steady influx 

175 of energy from the sun and re-radiation to outer space. Indeed, all living organisms 

176 survive (in relatively steady states) by balancing input with output - of energy and 

177 matter of some form. Labeling these reservoirs as 'the medium' in which our system 

178 finds itself, we see the following scenario emerging. Though the medium-|-system 

179 combination is clearly evolving in time, the system may be small enough that it has 

180 arrived at a time- independent NESS. While the much larger, combined system may 

181 well be evolving according to a time-reversal symmetric dynamics, it is quite reasonable 

182 to assume that this symmetry is violated by the effective dynamics appropriate for our 

183 smaller system with its shorter associated time scales. In other words, when we sum 

184 over the degrees of freedom associated with the medium, the dynamics describing the 

185 remaining configurations C's of our system should in general violate detailed balance. 

186 In general, it is impossible to derive such effective dynamics for systems at the 

187 mesoscopic or macroscopic scales from well-known interactions at the microscopic, 

188 atomic level. There are proposals to derive them from variational principles, based 

189 on postulating some quantity to be extremized during the evolution, in the spirit of 

190 least action in classical mechanics. The most widely known is probably 'maximum 

191 entropy production'. The major challenge is to identify the constraints appropriate 

192 for each NESM system at hand. None of these approaches has achieved the same level 

193 of acceptance as the maximum entropy principle in EQSM (where the constraints are 

194 well established, e.g., total energy, volume, particle number, etc.). In particular, the 

195 NESS in the uniformly driven lattice gas is known to differ from the state predicted 

196 by this principle [36]. Readers interested in these approaches may study a variety of 

197 books and reviews which appeared over the years [37l|38l|39l|40l|4ll|42l|43l|44]. How 

198 an effective dynamics (i.e., a set of rules of evolution) arise is not the purpose of our 

199 review. Instead, our goal is to explore the nature of stationary states, starting from 

200 a given dynamics that violates detailed balance. Specifically, in modeling biological 

201 transport, the main theme of the applications section, it is reasonable to postulate a 

202 set of transition rates for the system of interest. 
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203 2.1. Master equation and other approaches to statistical mechanics 

204 Since probabilities are central to statistical mechanics, our starting point for discussing 

205 NESM is P{C,t) and its evolution. Since much of our review will be devoted to 

206 models well-suited for computer simulations, let us restrict ourselves to discrete steps 

207 (r = 0, 1, ... of time 5t) as well as a discrete, finite set of C's (Ci, C2, .., C^f). Also, 

208 since we will be concerned with time reversal, we will assume, for simplicity, that 

209 all variables are even under this operation (i.e., no momenta- like variables which 

210 change sign under time reversal). To simplify notation further, let us write Pi (r), 

211 interchangeably with P {Ci,T6t). Being conserved (X^i ^« (''') — 1 ''')' ^ must 

212 obey a continuity equation, i.e., the vanishing of the time rate of change of a conserved 

213 density plus the divergence of the associated current density. Clearly, the associated 

214 currents here are probability currents. Since we restrict our attention to a discrete 

215 configuration space, each of these currents can be written as the flow from Cj to C^, 

216 i.e., Kf (t) or K {Cj — )■ Ci^rSt). As a net current, iiTj (r) is by definition —Kl (r), 

217 while its 'divergence' associated with any Ci is just 'Y^- (r). In general, K is a. new 

218 variable and how it evolves must be specified. For example, in quantum mechanics, P 

219 is encoded in the amplitude of the wavefunction -0 only, while K contains information 

220 of the phase in -0 as well (e.g., iiT cx -0* V for a single non-relativistic particle). In 

221 this review, as well as the models used in essentially all simulation studies, we follow 

222 a much simpler route: the master equation or the Markov chain. Here, K is assumed 

223 to be proportional to P, so that Pi (r -I- 1) depends only on linear combinations of 

224 the probabilities of the previous step, P, (r). In a further simplification, we focus on 

225 time-homogeneous Markov chains, in which the matrix relating X to P is constant in 

226 time. Thus, we write Pi (r + 1) = wjPi (t), where are known as the transition 

227 rates (from Cj to Ci). 

228 As emphasized above, we will assume that these rates are given quantities, as in a 

229 mathematical model system like TASEP or in phenomenologically-motivated models 

230 for biological systems. Probability conservation imposes the constraint J2i "^i — ^ for 

231 all j, of course. A convenient way to incorporate this constraint is to write the master 

232 equation in terms of the changes 

AP,(r) =P,(r+l)-P,(T) 

233 This equation can be written as 

A/^(r) = ^L^P,(r) (2) 
i 

235 where 

237 is a matrix (denoted by L; sometimes referred to as the Liouvillian) that plays much the 

238 same role as the Hamiltonian in quantum mechanics. Since all transition rates are non- 
239 negative, is a stochastic matrix and many properties of the evolution of our system 

240 follow from the Perron- Frobenius theorem [l^. In particular, J^i^l — for all j 

241 (probability conservation), so that at least one of the eigenvalues must vanish. Indeed, 
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242 we recognize the stationary distribution, P*, as the associated right eigenvector, since 

243 LP* — imphes P* (r + 1) = P* (r). Also, this P* is unique, provided the dynamics 

244 is ergodic, i.e., every Ci can be reached from any Cj via the w's. Further, the real 

245 parts of all other eigenvalues must be negative, so that the system must decay into 

246 P* eventually. 

247 Since the right-hand side of equation ([T]) is already cast in the form of the 

248 divergence of a current, we identify 

KI{t)^wIP,{t)^w)P,(t) (4) 

250 as the (net) probability current from Cj to d. Note that the antisymmetry Kl = —Kj 

251 is manifest here. When a system settles into a stationary state, these time-independent 

252 currents are given simply by 

Kr^wiP*-w)P: (5) 

254 As we will present in the next subsection, a reasonable way to distinguish EQSS from 

255 NESS is whether the K* vanish or not. Before embarking on that topic, let us briefly 

256 remark on two other common approaches to time dependent statistical mechanics. 

257 More detailed presentations of these and related topics are beyond the scope of this 

258 review, but can be found in many books [IHl 271 SHI HH] ■ 

259 Arguably the most intuitive approach to a stochastic process is the Langevin 
equation. Originating with the explanation of Brownian motion |5Qj by Einstein and 
Smoluchowski [Ml HI] , this equation consists of adding a random drive to an otherwise 
deterministic evolution. The deterministic evolution describes a single trajectory 
through configuration space: C (r) (starting with C (0) = Co), governed by say, an 
equation of the form AC (t) — J^[C{t)]. In a Langevin approach, T will contain 
both a deterministic part and a noisy component. Of course, a trajectory (or history, 

266 or realization) will depend on the specific noise force appearing in that run. Many 

267 trajectories are therefore generated, each depending on a particular realization of 
the noise and the associated probability. Although each trajectory can be easily 
understood, the fact that many of them are possible means the system at time r 

270 can be found at a collection of C's. In this sense, the evolution is best described by 

271 a probability distribution, P(C,t), which is controlled by both the deterministic and 

272 the noisy components in T. Historically, such considerations were first provided for a 

273 classical point particle, where AC = would be Newton's equation, dfx {t) = F/m, 

274 with continuous time and configuration variables. How the deterministic and noisy 

275 parts of F are connected to each other for the Brownian particle is the celebrated 

276 Einstein- Smoluchowski relation. Of course, P(C,t) becomes P{x,t) in this context, 

277 while the Langevin approach can be reformulated as a PDE for P {x, t) 

dtPix,t)= D^p{x)P{x,t)- —Va,{x)P{x,t) (6) 

279 Here, Da/s and Va are the diffusion tensor and the drift vector, respectively, and are 

280 related to the noisy and deterministic components of the drive[^. This PDE is referred 

281 to as the Fokker-Planck equation, although it was first introduced for the velocity 

282 distribution of a particle [15]. 

S Note that Dapix) here derives from the rate of change of the variance of the distribution and is 
different from the diffusion coefficient used to define Fick's law. Here, both spatial derivatives operate 
on Dc,i3{x). 



260 
261 
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283 An experienced reader will notice that the master equation ([T} for P (C,t) and the 

284 Fokker-Planck equation ([5]) for P {x, t) are both linear in P and first order in time, 

285 but that the right-hand sides are quite different. Let us comment briefly on their 

286 similarities and differences. Despite the simpler appearance, equation ^ is the more 

287 general case, apart from the complications associated with discrete vs. continuous 

288 variables. Thus, let us facilitate the comparison by considering a discrete version of 

289 equation (|6]), i.e., t — > rSt and x — >■ (6x, so that P{x,t) — >• P(<^,t). In this light, 

290 it is clear that the derivatives on the right correspond to various linear combinations 

291 of P{Ca ± Ij C/3 i 1; ''')• In other words, only a handful of the 'nearest configurations' 

292 are involved in the evolution of P(C,r). By contrast, the range of wf, as written in 

293 equation ([T]), is not restricted. 

294 Let us illustrate by a specific example. Consider a system with N Ising spins (with 

295 any interactions between them) evolving according to random sequential Glauber spin- 

296 flip dynamics 53J. In a time step, a random spin is chosen and flipped with some 

297 probability. Now, as noted earlier, the configuration space is the set of vertices of 

298 an N dimensional cube. Therefore, the only transitions allowed are moves along an 

299 edge of the cube, so that the range of is 'short'. In this sense, field theoretic 

300 formulations of the Ising system evolving with Glauber dynamics are possible, taking 

301 advantage of Fokker-Planck like equations, cast in terms of path integrals. On the 

302 other hand, consider updating according to a cluster algorithm, e.g., Swendsen-Wang 

303 [54] ■ in which a large cluster of spins (say, M) are flipped in a single step. Such a 

304 move clearly corresponds to crossing the body diagonal of an M dimensional cube. 

305 Since M is conceivably as large as N, there is no limit to the range of this set of w's. 

306 It is hardly surprising that fleld theoretic approach for such systems are yet to be 

307 formulated, as they would be considerably more complex. 

308 2.2. N on- equilibrium vs. equilibrium stationary states, persistent probability currents 

309 Following the footsteps of Boltzmann and Gibbs, we study statistical mechanics of 

310 systems in thermal equilibrium by focusing on time independent distributions such 

311 as P* (C) cx 1 or exp [~PT-L (C)] (where T-L is the total energy associated with C and 

312 /3 is the inverse temperature). Apart from a few model systems, it is not possible to 

313 compute, analytically and in general, averages of observable quantities, i.e., 

{0)^Y.'^{Cj)P*{C,). (7) 

315 Instead, remarkable progress over the last fifty years was achieved through computer 

316 simulations, in which a small subset of {C} is generated - with the appropriate 

317 (relative) weights - and used for computing the desired averages. This approach is an 

318 advanced art 55], far beyond the scope (or purpose) of this review. Here, only some 

319 key points will be mentioned and exploited - for highlighting the contrast between the 

320 stationary distributions of Boltzmann-Gibbs and those in NESS. 

321 In a classic paper [56], Metropolis, et.al. introduced an algorithm to generate 

322 a set of configurations with relative Boltzmann weights. This process also simulates 

323 a dynamical evolution of the system, in precisely the sense of the master equation. 

324 Starting from some initial C (0) = Cq, a new one, Cfc, is generated (by some well defined 

325 procedure) and accepted with probability Thus, C (1) is Ck or Cq. with relative 

326 probability w^/ (l — w^) . After some transient period, the system is expected to settle 

327 into a stationary state, i.e., the frequencies of Ci occurring in the run are proportional 
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to a time independent P* (Cj). To implement this Monte-Carlo method, a set of 
transition rates, w*, must be fixed. Further, if the desired outcome is P* cx e^*^^, 
then cannot be arbitrary. A sufficient (but not necessary) condition is referred to, 
especially in the simulations community |55| . as 'detailed balance': 

wlP* (C) = w^P* [Cu) . (8) 

In other words, it suffices to constrain the ratio w\/w^ to be exp[— /JAH], where 
A H = H (Cfc) — H (Ci) is just the difference between the configurational energies. A 
common and simple choice is wl. — min {l, e~^^^} . 

Of course, constraining the ratios still leaves us with many possibilities. To narrow 
the choices, it is reasonable to regard a particular set of w's as the simulation of a 
physical dynamicfli|. In that case, other considerations will guide our choices. For 
example, the Lenz-Ising system is used to model spins in ferromagnetism j57j as well 
as occupations in binary alloys |58[ 159] . In the former, individual spins can be flipped 
and it is quite appropriate to exploit Glauber |53] dynamics, in which the w's connect 
C's that differ by only one spin. In the latter however, a Zn atom, say, cannot be 
changed into a Cu atom, so that exchanging a neighboring pair of different 'spins' - 
Kawasaki |601 161] dynamics - is more appropriate. In terms of the N dimensional 
cube representation of {C}, these w's connect two vertices along an edge (Glabuer) or 
across the diagonal of a square face or plaquette (Kawasaki) . Both dynamics involve 
w's that only connect C's with one or two different spins. The idea is that, in a short 
6t, exchanging energy with the thermal reservoir can randomly affect only one or two 
spins. Also, in this sense, we can regard the w's as how the system is coupled to 
the surrounding medium. Clearly, AH measures the energy exchanged between the 
two. Another important quantity is entropy production, whether associated with the 
system or the medium, in which the w's will play a crucial role. 

It is significant that regardless of the details of the associated dynamics, a set of 
w's that obey detailed balance ([S]) necessarily leads the system to a state in which all 
stationary currents vanish. This follows trivially from the definition ((Sj. By contrast, 
transition rates that violate detailed balance necessarily lead to some non-vanishing 
stationary currents. To appreciate this statement, let us provide a better criterion, 
due to Kolmogorov [62], for rates that respect/ violate detailed balance. In particular, 
while equation gives the wrong impression that detailed balance is defined with 
respect to a given H, the Kolmogorov criterion for detailed balance is applicable to 
all Markov processes, whether an underlying Hamiltonian exists or not. 

Consider a closed loop in configuration space, e.g., £ = —> Cj —>■ 
Cn — ^ Ci. Define the product of the associated rates in the 'forward' direction by 
Il[C] = Wjwl. ...w"^ and also for the 'reverse' direction: 11 [£rei>] = •■•^n- The 

set of rates are said to satisfy detailed balance if and only if 

n [£] = n [Crev] (9) 

for all loops. If this criterion is satisfied, then we can show that a (single valued) 
functional in configuration space can be constructed simply from the set of ratios 
w^/w^, and that it is proportional to P* (C). If this criterion is violated for certain 
loops, these will be referred to here as 'irreversible rate loops' (IRLs). Despite the lack 

II In this light, the sequence of configurations generated {Cj-^yCj^, . . . ,Cj^) can be regarded as a 
history, or trajectory, of the system. By collecting many (M) such sequences, | C" j- ,a = 1, ...,M, 

time dependent averages (O)^ = ^ (Cj) Pj (t) are simulated by ^ i^j' 
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371 of detailed balance, P* (C) exists and can still be constructed from the w's, though 

372 much more effort is required. Established some time ago [31|SS1[B1], this approach to P* 

373 is similar to Kirchhoff's for electric networks |65| . More importantly, this construction 

374 provides the framework for showing that, in the stationary state, some K*^s must be 

375 non-trivial [66l [7] . Since the divergence of these currents must vanish, they must form 

376 current loops. As time-independent current loops, they remind us of magnetostatics. 

377 The distinction between this scenario and the case with detailed balance w's is clear: 

378 The latter resembles electrostatics. In this light, it is reasonable to label a stationary 

379 state as an equilibrium one if and only if all its (probability) currents vanish, associated 

380 with a set of w's with no IRLs. Similarly, a non-equilibrium steady state - NESS - 

381 would be associated with non-trivial current loops, generated by detailed balance- 

382 violating rates with IRLs [SSI H] ■ Our proposal is that all stationary states should be 

383 characterized by the pair {P*,if*}. In this scheme, 'equilibrium states' correspond 

384 to the subset {P*,0}, associated with a dynamics that respect detailed balance and 

385 time reversal. 

386 The presence of current loops and IRLs raises a natural and intriguing question: 

387 Is there a intuitively accessible and simple relationship between them? Unfortunately, 

388 the answer remains elusive so far. Venturing further, it is tempting to speculate on the 

389 existence of a gauge theory, along the lines of that in electromagnetism, for NESM. If 

390 such a theory can be formulated, its consequences may be far-reaching. 

391 Time- independent probability current loops also carry physical information about 

392 a NESS. Referring readers to a recent article [7] for the details, we provide brief 

393 summaries here for a few key points. 

394 (i) In particular, it is shown how the if's can be used to compute currents associated 

395 with physical quantities, such as energy or matter. In addition, we have emphasized 

396 that a signature of NESS is the existence of a steady flux (of, e.g., energy) through the 

397 system. All aspects of such through-flux, such as averages and correlation, can also 

398 be computed with the i^'s. 

399 (ii) Following Schnakenberg [63], we may define the total entropy production Stot as 

400 a quantity associated with the rates { } . This Stot can be written as the sum of two 

401 terms, Sgyg-I- Smcd- The first is associated with entropy production within our system 

402 (recognizable as the derivative of the Gibbs' entropy, — Pi InP^, in the continuous 

403 time limit): 

S.ys(r)^i5]A^^"(r)ln^. (10) 

405 It is straightforward to show that for a NESS with K* ^ 0, Sgyg vanishes as expected. 

406 However, a second contribution to entropy production is associated with the medium: 

S„ed(r)^i5]X/(r)ln4, (11) 

408 where the positivity of ^ AT*-* ln(w^ / ) has been demonstrated. This result is entirely 

409 consistent with our description of a NESS, namely, a system coupled to surroundings 

410 which continue to evolve and generate entropy. 

411 (iii) The following inverse question for NESS is also interesting. As we noted, given a 

412 Boltzmann distribution, a well known route to generate it is to use a dynamics which 
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413 obey detailed balance ([5]). If we accept that a NESS is characterised not only by the 

414 stationary distribution, but by the pair {P*,K*}, then the generalized condition is 

415 wl.P* = wfP^ + K", or more explicitly, 

wlP*iC,)^w^P*iCk)+K*{C^^Ck). (12) 

417 It is possible to phrase this condition more elegantly (perhaps offering a little insight) 

418 by performing a well- known similarity transform on wl.: Define the matrix U, the 

419 elements of which are 

ui^ip^r'^'wiip*)'^'. 

421 The advantage of this form of 'coding' the dynamics is that U is symmetric if rates 

422 obey detailed balance. Meanwhile, since K* is a current, we can exploit the analog 

423 J — pv to define the 'velocity matrix' V, the elements of which are 

^ {p^r'/' Kr{p*r'^\ 

425 associated with the flow from Ci to Ck- Our generalized condition (I12[) now reads 

426 simply: The antisymmetric part of U is V/2. Similar ideas have also been pursued 

427 recently [67 . 

428 To summarize, if a dynamics is to lead a system to a desired {P*,K*}, then 

429 the associated antisymmetric part of M is completely fixed by K*. By contrast, its 

430 symmetric part is still unconstrained, corresponding to dynamics that takes us to the 

431 same {P*,K*}. 

432 (iv) As long as K* 7^ for a transition, we can focus on the direction with positive 

433 current (say, K^^ > 0) and choose the maximally asymmetric dynamics, namely, 

434 = and = K^'^/P* . (Understandably, such choices are impossible for systems in 

435 thermal equilibrium, except for T = cases.) Whether systems with such apparently 

436 unique dynamics carry additional significance remains to be explored. Certainly, 

437 TASEP - the paradigmatic model of NESS, to be presented next - belongs in this 

438 class. Before embarking on the next section, let us briefiy comment on some typical 

439 features of NESS which are counter-intuitive, based on our notions of EQSM. 

440 2. 3. Beyond expectations of equilibrium statistical mechanics 

441 Equilibrium statistical mechanics has allowed us to develop physical intuition that can 

442 be valuable guides when we are faced with new problems in unfamiliar settings. A 

443 good example is energy-entropy competition, which tends to serve us well when we 

444 encounter novel phase transitions: The former/latter 'wins' for systems at low/high 

445 temperatures, so that it displays order/disorder phenomena. Another example is 

446 'positive response': To ensure thermodynamic stability, we expect the system to 

447 respond in a certain manner (positively) when its control parameters are changed. 

448 Thus, it is reasonable to expect, e.g., positive specific heat and compressibility for 

449 systems in thermal equilibrium. A final example is long-range correlations, which are 

450 generically absent in equilibrium systems with short-range interactions and dynamics. 

451 There are exceptions, of course, such as in critical phenomena associated with second 

452 order phase transitions. When we encounter systems in NESS however, we should be 

453 aware that such physical intuition often leads us astray. At present, we are not aware of 

454 another set of overarching principles which are generally applicable for NESS. Instead, 

455 in the following, we will provide two specific circumstances in which our expectations 

456 fail. 
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457 Negative response. In order for a system to be in thermal equilibrium, it must be 

458 stable against small perturbations. Otherwise, fluctuations will drive it into another 

459 state. Such consistent behaviors of a system may be labeled as 'positive response'. 

460 Related to the positivity of certain second derivatives of the free energy, elementary 

461 examples include positive specific heat and compressibility. By contrast, a surprisingly 

462 common hallmark of NESS is 'negative response.' For example, imagine a room in 

463 which the internal energy decreases when the thermostat is turned up! One of the first 

464 systems in NESS where this type of surprising behavior surfaced is the driven Ising 

465 lattice gas Referring the reader to, for example, reference [55] for details, the 

466 key ingredients are the following. An ordinary Ising system (with nearest neighbor 

467 ferromagnetic interactions on a square lattice) is subjected to an external drive, and 

468 observed to undergo the phase transition at a temperature higher than that expected 

469 from Onsager's solution. Since the external drive tends to break bonds, its effect is 

470 similar to coupling the system to another energy reservoir with a higher temperature. 

471 Nevertheless, this NESS system displays more order than its equilibrium counterpart. 

472 In other words, despite the fact that more energy appears to be 'pumped into' the 

473 system, the internal energy decreases. A more direct manifestation of this form of 

474 'negative response' has been observed in the two-temperature Ising lattice gas, in 

475 which particle hops in the x or y direction are updated by Metropolis rates appropriate 

476 to exchanging energy with a thermal reservoir set at temperature or Ty. Changing 

477 Tx, with Ty held fixed, the average internal energy, U (i.e., (H)), is found to vary with 

478 dU/ dTy < [69] ! Such surprising negative response is so generic that it can be found 

479 in exceedingly simple, exactly solvable cases [70] . 

480 Of course, we should caution the reader that 'negative response' may be simply 

481 a misnomer, poor semantics, or careless interpretation of an observed phenomenon. 

482 After all, fluctuations of observables in a stationary state must be positive and if 

483 the appropriate conjugate variable is used, then the response to that variable will 

484 again be positive. In particular, for any observable O, we can always define the 

485 cumulant generating function S (w) = ln(e'^'^) and its derivative X (uj) = E' {uj). 

486 Of course, the average (O) is S'(0), while (i^) is, in general, the average of 

487 O in the modified distribution P* (C) cx e'^'^'^'^-'P* (C). Then, we are guaranteed 

488 'positive response:' dX/duj > 0. However, unlike internal energy and temperature for 

489 systems in thermal equilibrium, simple physical interpretations of these mathematical 

490 manipulations for NESS's are yet to be established. Clearly, this issue is related 

491 to the fiuctuation-dissipation theorem in EQSM. General results valid for NESM 

492 have been derived during the last two decades; we refer the reader to the seminal 

493 articles [71] [72] [73] [27] [28] . The generalization of the fluctuation-dissipation theorem 

494 to NESM is a major topic [741 134] and lies outside our scope. Here, let us remark that 

495 the foundations of this theorem lies in the time reversal properties of the underlying 

496 dynamics [75] [76], which control the nature of the fluctuations of the random variables. 

497 To characterize these fluctuations quantitatively, large-deviation functions (LDF) have 

498 been introduced. They play a crucial role in NESM, akin to that of thermodynamic 

499 potentials in EQSM [77] [78] . Valid for systems far from equilibrium, the fluctuation 

500 theorem can be stated as a symmetry property of the LDF (see the next section for 

501 a explicit example in the case of TASEP). Near an equilibrium state this theorem 

502 implies the fluctuation-dissipation relation, previously derived from linear response 

503 theory [7n][H0]- A related set of significant results is the non-equilibrium work relations 

504 [27] [3I] [29] [2 81 [301 [351 [331 134l |3T] ■ also a topic worthv of its own re vie w ( see, e . g. . [8n[82] 

505 and references therein) . Here, in the context of the exact solution of TASEP (section 
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506 [3]), another result along this theme - the macroscopic fluctuation theory developed by 

507 Jona-Lasinio et al. [8 3) - plays an important role. 

508 Generic long-range correlations. For systems with short-range interactions in thermal 

509 equilibrium, the static (i.e., equal time) correlations of the fluctuations are generically 

510 short-ranged. This is true even if the dynamics were to contain long ranged 

511 components (e.g., the 'artificial' Swendsen-Wang [54] dynamics): As long as it 

512 obeys detailed balance, we are assured of P* cx e~^^. Of course, time-dependent 

513 correlations are not similarly constrained. An excellent example is diffusive dynamics 

514 obeying certain conservation laws (e.g., hydrodynamics or Kawasaki [60[I61| dynamics 

515 modeling Cu-Zn exchange), where long time tails (power law decays) are well 

516 known phenomena: The autocorrelation function, G {r = 0,t), decays as t~'^^'^ in d 

517 dimensions, even for a single particle. As pointed out by Grinstein j84j . a simple 

518 scaling argument, along with r ~ t^^^ , would lead to generic long-range correlations, 

519 G {r,t — 0) Ar^'"^^^ (i.e., r^'^ in the case of random walkers subjected to short- 

520 range interactions, where z = 2). Yet, G {r,t = 0) generally decays as an exponential 

521 in equilibrium systems! These seemingly contradictory statements can be reconciled 

522 when the amplitude A is examined. In the equilibrium case, detailed balance (or the 

523 fluctuation-dissipation theorem) constrains A to vanish. For NESS, there is no such 

524 constraint, leaving us with long-range correlations generically. For further details of 

525 these considerations, the reader may consult references [551155] . 

526 To end the discussion on such correlations, we should caution the readers on a 

527 subtle point. Although we emphasized how long-range correlations can emerge from 

528 a short-range dynamics that violates detailed balance, the latter is not necessary. 

529 An excellent example is a driven diffusive system with three species - the ABC 

530 model [Sni [HI] ~ in one dimension. The system evolves with a short-range dynamics 

531 which generally violate detailed balance, and displays long-range order (as well as 

532 correlations). Remarkably, for a special set of parameters, detailed balance is restored 

533 and so, an exact P* was easily found. When interpreted as a Boltzmann factor (i.e., 

534 P* cx e~^'^), the Hamiltonian contains inter-particle interactions which range over 

535 the entire lattice! Despite having an Ji with long-range interactions, it is possible 

536 to construct a short-range dynamics (e.g., w's that involve only nearest neighbor 

537 particle exchanges) that will lead the system to an EQSS: {P* cx e"^^, K* = O}. To 

538 appreciate such a counter-intuitive situation, consider AJi — 'H{Ci)— 'H{Cj) for a 

539 pair of configurations that differ in some local variables. The presence of long-range 

540 terms in T-L typically induce similar terms in AH, leading to a long-range dynamics. 

541 If such terms conspire to cancel, then AH becomes short-ranged and it is simple to 

542 choose u;'s with no long range components. We believe it is important to investigate 

543 whether such examples belong to a class of mathematical curiosities or form the basis 

544 for a wide variety of physical/biological systems. With these two examples, we hope to 

545 have conveyed an important lesson we learned from our limited explorations of NESS: 

546 Expect the unexpected, when confronted with a novel system evolving according to 

547 detailed balance violating dynamics. 

548 In this section, we attempted to provide a bird's eye view of the many facets 

549 of non-equilibrium statistical mechanics, and then to focus on a particular aspect: 

550 stationary states associated with dynamics that violate detailed balance. We 

551 emphasized the importance of this class of problems and pointed out significant 

552 features of NESS that run counter to the physical intuition learned from equilibrium 

553 statistical mechanics. In the next two sections, we will turn our attention to specific 

554 systems. As common in all theoretical studies, there are two typically diverging goals 
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555 associated with the models we pursue. One goal is to account for as many features 

556 found in nature as reaUstically possible. The other is to consider models with just one 

557 or two features so that they are simple enough to be solved analytically and exactly. 

558 These goals obviously diverge since more realistic models are typically mathematically 

559 intractable while exactly solvable models generally lack essential aspects of physical 

560 systems (or those in chemistry, biology, finance, sociology, etc.). Nevertheless, we 

561 believe is it worthwhile to devote our attention to both of them, albeit separately. In 

562 this spirit, we will next present a simple model: the exclusion process, with the TASEP 

563 being an extreme limit. Arguably the simplest possible model with a nontrivial NESS, 

564 the TASEP is not only amenable to exact and approximate solution strategies, but it 

565 also has shed considerable light on problems in the real world. In section ([4]), we turn 

566 to a number of generalizations of this model, each taking into account new physical 

567 features required for modeling various aspects of transport in biological systems. 

568 3. A paradigmatic model: The asymmetric simple exclusion process 

569 (ASEP) 

570 Building a simple representation for complex phenomena is a common procedure in 

571 physics, leading to the emergence of paradigmatic models: the harmonic oscillator, the 

572 random walker, the Ising magnet. All these beautiful models often display wonderful 

573 mathematical structures [551 In the field of NESM, and for investigating the role 

574 of detailed balance violating dynamics in particular, the asymmetric simple exclusion 

575 process (ASEP) is reaching the status of such a paradigm. A model with possibly 

576 the simplest of rules, it nevertheless displays a rich variety of NESS's. Further, it is 

577 sufficiently simple to allow us to exploit rigorous mathematical methods and, over the 

578 last two decades, to arrive at a number of exact results. In this way, valuable insights 

579 into the intricacies of NESM have been garnered. In this section, we delve into some 

580 details of two such methods, in the hope that readers unfamiliar with these techniques 

581 can join in these efforts. Of course, as we consider models which are more suited 

582 for physical/biological systems, we will encounter more complex ingredients than in 

583 ASEPs. As a result, these models are not exactly solvable at present. In these cases, 

584 approximations are necessary for further progress. One successful scheme is the mean 

585 field approximation, leading to hydrodynamics equations (PDE's) for density fields. 

586 Since this strategy is the only one to offer some quantitative understanding of the 

587 more realistic/complex models, we will devote the last subsection 13.111 here to this 

588 approach. 

589 In the previous section, we presented the master equation ([1]) in discrete time, 

590 which is clearly the most appropriate description for computer simulations. On the 

591 other hand, for analytic studies, it is often far easier to use continuous variables (or 

592 infinite systems, in a similar vein). Thus, all the discussions in this section will be 

593 based on continuous time, t. As discussed in the context of the Fokker-Planck equation 

594 ((6]), we introduced this connection: t = rSt. Here, let us be explicit and write the 

595 continuous version of equation ^ as 



597 where the matrix on the right, M, is just h/St. Taking the limit (5< — > in 

598 the formal solution, P (r) = [I + L]^P(0) (I being the identity), leads then to 

599 P{t) =exp[Mt]P{0). 



596 




(13) 
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600 3.1. Motivation and definition of ASEP and TASEP 

601 The ASEP is a many-body dynamical system, consisting of particles located on a 

602 discrete lattice that evolves in continuous time. The particles hop randomly from a 

603 site on a lattice to one of its immediate neighbors, provided the target site is empty. 

604 Physically, this constraint mimics short-range interactions amongst particles. In order 

605 to drive this lattice gas out of equilibrium, non- vanishing currents must be established 

606 in the system. This can be achieved by various means: by starting from non-uniform 

607 initial conditions, by coupling the system to external reservoirs that drive currents 

608 [90] through the system (transport of particles, energy, heat), or by introducing 

609 some intrinsic bias in the dynamics that favors motion in a privileged direction. 

610 Then, each particle is an asymmetric random walker that drifts steadily along the 

611 direction of an external driving force. Due to its simplicity, this model has appeared 

612 in different contexts. As noted in the Introduction, it was first proposed by Gibbs, et 

613 al. in 1968 [131 HI] as a prototype to describe the dynamics of ribosomes translating 

614 along an mRNA. In the mathematical literature, Brownian processes with hard-core 

615 interactions were defined in 1970 by Spitzer l9| who coined the term 'exclusion process' 

616 (see also [9Tl|92j[TTl[T0]). In addition to motivation from issues in molecular biology - 

617 the main focus of section lU the ASEP has also been used to describe transport in low- 

618 dimensional systems with strong geometrical constraints 1931 such as macromolecules 

619 transiting through anisotropic conductors, or quantum dots, where electrons hop to 

620 vacant locations and repel each other via Coulomb interaction [94], while many of 

621 its variants are ubiquitous in modeling traffic flow [95] [15]. More generally, the 

622 ASEP belongs to the class of driven diffusive systems defined by Katz, Lebowitz and 

623 Spohn in 1984 [Ml IMl [111211 [S3] ■ We emphasize that the ASEP is defined through 

624 dynamical rules, while no energy is associated with a microscopic configuration. More 

625 generally, the kinetic point of view seems to be a promising and fruitful approach to 

626 non-equilibrium systems (see e.g., [97]). 



627 To summarize, the ASEP is a minimal model to study non-equilibrium behavior. 

628 It is simple enough to allow analytical studies, however it contains the necessary 

629 ingredients for the emergence of a non-trivial phenomenology: 

630 • Asymmetric: The external driving breaks detailed balance and creates a 

631 stationary current. The system settles into a NESS. 

632 • Exclusion: The hard-core interaction implies that there is at most one particle 

633 per site. The ASEP is a genuine many-body problem, with arguably the simplest 

634 of all interactions. 

635 • Process: With no underlying Hamiltonian, the dynamics is stochastic and 

636 Markovian. 

637 Having a general picture of an ASEP, let us turn to a complete definition of the 



638 model. Focusing on only exactly solvable cases, we restrict ourselves here to a one 

639 dimensional lattice, with sites labeled by i = 1,...,L (here, we will use i to label a 

640 site rather than a configuration). The stochastic evolution rules are the following: A 

641 particle located at site i in the bulk of the system jumps, in the interval [t, t + dt], with 

642 probability pdt to site i + 1 and with probability qdt to site i — 1, provided the target 

643 site is empty {exclusion rule). The rates p and q are the parameters of our system. By 

644 rescaling time, p is often set to unity, while q is left as a genuine control parameter. 

645 In the totally asymmetric exclusion process (TASEP), the jumps are totally biased in 

646 one direction (e.g., g = 0). On the other hand, the symmetric exclusion process (SEP) 
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647 corresponds to the choice p ~ q. The physics and the phenomenology of the ASEP is 

648 extremely sensitive to the boundary conditions. We wih mainly discuss three types of 

649 boundary conditions: 

650 Periodic boundary conditions: 

651 The exclusion process is defined on a ring, so that sites i and L + i are identical. 

652 The system is filled with N < L particles (figure [TJa)). Of course, the dynamics 

653 conserves N. 

654 Infinite lattice: 

655 Here, the boundaries are sent to ±oo. Boundary conditions are here of a different 

656 kind. This system always remains sensitive to the initial conditions. Therefore, an 

657 initial configuration, or a statistical set of configurations, must be carefully specified. 

658 Figure [Ijb) is an illustration of this system. 

659 Open boundaries: 

660 Here, the boundary sites i — 1,L play a special role, as they are coupled to 

661 particles in reservoirs just outside the system. Thus, if site 1 is empty, a particle can 

662 enter (from the 'left reservoir') with rate a. If it is occupied, the particle can exit 

663 (into this reservoir) with rate 7. Similarly, the interactions with the 'right reservoir' 

664 are: If site L is empty/occupied, a particle can enter/exit the system with rate S//3 

665 respectively. These rates can be regarded as the coupling of our finite system with 

666 infinite reservoirs set at different 'potentials'. Figure [Ijc) illustrates this system. A 

667 much simpler limiting case is the TASEP, with g = 7 = (5 = 0, in which particles 

668 injected from the left, hopping along the lattice only to the right, and finally exiting 

669 to the right. 

9 P 
• !•! !•! !•! I»l I > 

(b) 





Figure 1. (a) The asymmetric simple exclusion process on a ring of L sites filled 
with N particles (L = 24, A'^ = f2 here). The total number of configurations is 
just n = (^). (b) ASEP on an infinite lattice. The particles perform asymmetric 
random walks (p 7^ q) and interact through the exclusion constraint, (c) A 
schematic of a ASEP with open boundaries on a finite lattice with L = 10 sites. 

670 We emphasize that there are very many variants of the ASEP. The dynamical 

671 rules can be modified, especially in computer simulation studies using discrete-time 

672 updates (e.g., random sequential, parallel, or shuffle updates). The hopping rates 

673 can be modified to be either site- or particle-dependent, with motivations for such 

674 additions from biology provided below. In modeling vehicular traffic, the former is 

675 suitable, e.g., for including road work or traffic lights. The latter can account for the 

676 range of preferred driving velocities, while a system can also be regarded as one with 
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677 many 'species' of particles. In addition, these disorders can be dynamic or quenched 

678 [Ml ISH Enni [Ml [101 [ini lini ■ Further, the exclusion constraint can be relaxed, 

679 SO that fast cars are allowed to overtake slower ones, which are known as 'second class' 

680 or 'third class' particles [inSl[inZl[inSllinil[nni[ini[II2- Finally, the lattice geometry 

681 itself can be generalized to multi-lanes, higher dimensions, or complex networks. All 

682 these modifications drastically alter the collective behavior displayed in the system, as 

683 hundreds of researchers discovered during the last two decades. Though more relevant 

684 for applications, more realistic models cannot, in general, be solved exactly. As the 

685 primary focus of this section is exact solutions, we will focus only on the homogeneous 

686 systems presented above. These problems are amenable to sophisticated mathematical 

687 analysis thanks to a large variety of techniques: Bethe Ansatz, quadratic algebras, 

688 Young tableaux, combinatorics, orthogonal polynomials, random matrices, stochastic 

689 differential equations, determinantal representations, hydrodynamic limits etc. Each 

690 of these approaches is becoming a specific subfield that has its own links with other 

691 branches of theoretical physics and mathematics. Next, we will present some of these 

692 methods that have been developed for these three ideal cases. 

693 3.2. Mathematical setup and fundamental issues 

694 The evolution of the system is given by equation p3)) and controlled by the Markov 

695 operator M as follows. In order distinguish the two uses of z - label for configurations 

696 and for sites on our lattices, let us revert to using C for configurations. Then, this 

697 master equation reads 



699 As a reminder, the off-diagonal matrix elements of M represent the transition rates, 

700 which the diagonal part M{C,C) — — J2c'^c M{C',C') accounts for the exit rate from 

701 C. Thus, the sum of the elements in any given column vanishes, ensuring probability 

702 conservation. For a finite ASEP, {C} is finite and the Markov operator M is a 

703 matrix. For the infinite system, M is an operator and its precise definition needs 

704 more sophisticated mathematical tools than linear algebra, namely, functional analysis 

705 [92j Jilj . Unless stated otherwise, we will focus here on the technically simpler case of 

706 finite L and deduce results for the infinite system by taking L ~> oo limit formally. An 

707 important feature of the finite ASEP is ergodicity: Any configuration can evolve to any 

708 other one in a finite number of steps. This property ensures that the Perron- Frobenius 

709 theorem applies (see, e.g., [^1113) ). Thus, i? = is a non-degenerate eigenvalue of M, 

710 while all other eigenvalues have a strictly negative real part, Re(£') < 0. The physical 

711 interpretation of the spectrum of M is clear: The right eigenvector associated with 

712 the eigenvalue E = corresponds to the stationary state of the dynamics. Because all 

713 non-zero eigenvalues have strictly negative real parts, these eigenvectors correspond 

714 to relaxation modes, with — 1/Re(i?) being the relaxation time and lm{E) controlling 

715 the oscillations. 

716 We emphasize that the operator M encodes all statistical information of our 

717 system, so that any physical quantity can be traced ultimately to some property of M. 

718 We will list below generic mathematical and physical properties of our system that 

719 motivate the appropriate calculation strategies: 

720 (i) Once the dynamics is properly defined, the basic question is to determine the 

721 steady-state, P*, of the system, i.e., the right eigenvector of M with eigenvalue 



698 
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722 0. Given a configuration C, the value of tlie component P*{C) is known as tfie 

723 measure (or stationary weiglit) of C in tlie steady-state, i.e., it represents tfie 

724 frequency of occurrence of C in the stationary state. 

725 (ii) The vector P* plays the role of the Boltzmann factor in EQSM, so that all steady- 

726 state properties (e.g., all equal-time correlations) can be computed from it. Some 

727 important questions are: what is the mean occupation pi of a given site i? What 

728 does the most likely density profile, given by the function i — >■ pi, look like? Can 

729 one calculate density-density correlation functions between different sites? What 

730 is the probability of occurrence of a density profile that differs significantly from 

731 the most likely one? The last quantity is known as the large deviation of the 

732 density profile. 

733 (iii) The stationary state is a dynamical state in which the system constantly 

734 evolves from one micro-state to another. This microscopic evolution induces 

735 macroscopic fluctuations (the equivalent of the Gaussian Brownian fiuctuations at 

736 equilibrium). How can one characterize such fluctuations? Are they necessarily 

737 Gaussian? How are they related to the linear response of the system to small 

738 perturbations in the vicinity of the steady-state? These issues can be tackled 

739 by considering tagged-particle dynamics, anomalous diffusion, time-dependent 

740 perturbations of the dynamical rules |114) . 

741 (iv) As expected, the ASEP carries a finite, non-zero, steady-state particle current, 

742 J, which is clearly an important physical observable. The dependence of J on 

743 the external parameters of the system allows us to define different phases of the 

744 system. 

745 (v) The existence of a non-zero J in the stationary state implies the physical transport 

746 of an extensive number of particles, Q, through the lattice. The total number of 

747 particles Qt, transported up to time t is a random quantity. In the steady-state, 

748 the mean value of Qt is just Jt, while in the long time limit, the distribution 

749 of the random variable Qt/t — J represents exceptional fiuctuations (i.e., large 

750 deviations) of the mean-current. This LDF is an important observable that 

751 provides detailed properties of the transport through the system. While particle 

752 current is the most obvious quantity to study in an ASEP, similar questions can 

753 be asked of other currents and the transport of their associated quantities, such 

754 as mass, charge, energy, etc., in more realistic NESM models. 

755 (vi) The way a system relaxes to its stationary state is also an important characteristic. 

756 The typical relaxation time of the ASEP scales with the system size as L^, where 

757 z is the dynamical exponent. The value of z is related to the spectral gap of the 

758 Markov matrix M, i.e., to the largest — Re(i?). For a diffusive system, z = 2. 

759 For the ASEP with periodic boundary condition, an exact calculation leads to 

760 z = 3/2. More generally, the transitory state of the model can be probed using 

761 correlation functions at different times. 

762 (vii) The matrix M is generally a non-symmetric matrix and, therefore, its right 

763 eigenvectors differ from its left eigenvectors. For instance, a right eigenvector 

764 ipE corresponding to the eigenvalue E is defined as 

765 MV' = Eil^. (15) 

766 Because M is real, its eigenvalues (and eigenvectors) are either real numbers 

767 or complex conjugate pairs. However, M is generally asymmetric, so that its 
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768 left eigenvectors are different from its right eigenvectors. Powerful analytical 

769 techniques, such as the Bethe Ansatz, can be exploited to diagonalize M in some 

770 specific cases, providing us with crucial information on its spectrum. 

771 (viii) Solving the master equation analytically would allow us to calculate exactly 

772 the evolution of the system. A challenging goal is to determine the finite-time 

773 Green function (or transition probability) P(C,t|Co,0), the probability for the 

774 system to be in configuration C at time t, given that the initial configuration at 

775 time t = Q was Cq. Formally, it is just the C,Co element of the matrix exp [Mi] 

776 here. In principle, it allows us to calculate all the correlation functions of the 

777 system. However, explicit results for certain quantities will require not only the 

778 knowledge of the spectrum and eigenvectors of M, but also explicit evaluations of 

779 matrix elements associated with the observable of interest. 

780 The following sections are devoted to a short exposition of several analytical 



781 techniques that have been developed to answer some of these issues for the ASEP. 

782 3.3. Steady state properties of the ASEP 

783 Given a stochastic dynamical system, the first question naturally concerns the 

784 stationary measure. We will briefiy discuss the ASEP with periodic boundary 

785 conditions and the infinite line case. More details will be provided for the highly 

786 non-trivial case of the open ASEP. 

787 

788 Periodic boundary conditions: 

789 This is the simplest case, with the stationary measure being flat [55] • That the 

790 uniform measure is also stationary can be understood as follows. A given configuration 

791 consists of k clusters of particles separated by holes. A particle that leads a cluster can 

792 hop in the forward direction with rate p whereas a particle that ends a cluster can hop 

793 backwards with rate g; thus the total rate of leaving a configuration consisting of k 

794 clusters is k{p + q). Similarly, the total number of ancestors of this configuration (i.e., 

795 of configurations that can evolve into it) is also given by k{p + q). The fact that these 

796 two quantities are identical sufltices to show that the uniform probability is stationary. 

797 To obtain the precise value of P*, l/il, is also elementary. Since iV is a constant, we 

798 only need the total number of configurations for N particles on a ring with L sites, 

799 whichisjust ^2 = L!/[iV!(L-iV)!]. 

800 

801 Infinite lattice: 

802 For the exclusion process on an infinite line, the stationary measures are studied 

803 and classified in \92\ lllj. There are two one-parameter families of invariant measures. 

804 One family, denoted by Vp, is a product of local Bernoulli measures of constant density 

805 p: this means that each site is occupied with probability p. The other family is 

806 discrete and is concentrated on a countable subset of configurations. For the TASEP, 

807 this second family corresponds to blocking measures, which are point-mass measures 

808 concentrated on step- like configurations (i.e., configurations where all sites to the 

809 left/right of a given site are empty/occupied). 

810 

811 Open boundaries: 

812 Turning to the case of the ASEP on a finite lattice with open boundaries, we 

813 note that the only knowledge we have, without detailed analysis, is the existence of a 
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814 unique stationary measure (thanks to the Perron- Frobenius theorem), i.e., the vector 

815 P* with 2^ components. We emphasize again that finding P* is a non-trivial task 

816 because we have no a priori guiding principle at our disposal. With no underlying 

817 Hamiltonian and no temperature, no fundamental principles of EQSM are relevant 

818 here. The system is far from equilibrium with a non-trivial steady-state current that 

819 does not vanish for even large L. 

820 To simplify the discussion, we focus on the TASEP where particles enter from the 

821 left reservoir with rate a, hop only to the right and leave the system from the site L 

822 with rate /?. A configuration C can be represented by the binary string, (cti, . . . , ctl), 

823 of occupation variables: 1 if the site i is occupied and u,; = otherwise. Our 

824 goal is to determine P* (C), with which the steady-state current J can be expressed 

825 simply and exactly: 

J = a (1 - (ai)) = - a^)) = . . . = - a.+i)) = /3(ai). (16) 

827 Here, the brackets ( ) denote expectation values as defined in equation ([7]). Even if J 

828 were known somehow, this set of equations is not sufficient for fixing the (unknown) 

829 density profile pi = (ai) and the nearest neighbor (NN) correlations (aiai^i). Typical 

830 in a truly many-body problem, there is a hierarchy |1151lll61lll71IT7] of equations that 

831 couple fc-site and (fc + l)-site correlations. A very important approximation, which 

832 often provides valuable insights, is the mean-field assumption in which the hierarchy 

833 is simply truncated at a given level. If the correlations beyond this level are small, this 

834 approximation can be quite good. Applying this technique here consists of replacing 

835 two-sites correlations by products of single-site averages: 

836 (c^iO-j) = PiPj. (17) 

837 Thus, equation (fT6| becomes 

838 J ~ a (1 - Pi) = pi(l - P2) = • • ■ = /5i(l - /Oj+i) = /3/3L- (18) 

839 and we arrive at a closed recursion relation between pi^i and p^, namely pi+i = 

840 1 — J/pi. Of course, J is an unknown, to be fixed as follows. Starting with pi = 1 — J/a, 

841 the recursion eventually leads us to pL as a rational function of J (and a). Setting 

842 this to J//3 gives a polynomial equation for J. Solving for J, we obtain the desired 

843 dependence of the steady-state current on the control parameters: J {a, (3, L). This 

844 approach to an approximate solution was known to Gibbs, et al. \13\ 114] (within the 

845 context of a more general version of TASEP) and explored further in [118j recently. 

846 Analyzing J {a, [3, L) gives rise to the phase diagram of the TASEP (see figure 

847 [5]) . For studying these aspects of the TASEP, the mean-field method provides us with 

848 a reasonably good approximation. Indeed, the correct phase diagram of the model 

849 was obtained in '90' \^ through physical reasoning by using such mean-field arguments 

850 along with the hydrodynamic limit. Since this strategy is quite effective and more 

851 widely applicable, we will briefiy discuss how it is applied to ASEP in section 13.111 

852 below. Despite many effective MFT-based strategies, exact solutions to ASEP are 

853 desirable, especially for an in-depth analysis. In particular, MET cannot account 

854 for correlations, fiuctuations, or rare events. In fact, the mean-field density profile 

855 (from solving the recursion relation ()18|) numerically) does not agree with the exact 

856 profile (from evaluating the expression (|23|) below). Of course, it is rare that we 

857 are able to calculate the stationary measure for a non-equilibrium interacting model. 

858 Not surprisingly, the exact steady-state solution of the TASEP jll8[ 11201 1121] was 

859 considered a major breakthrough. 

^ See also |10| and |119| for a pedagogical example. 
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860 3-4-. The matrix representation method for the open TASEP 

861 The exact calculation of the stationary measure for the TASEP with open boundaries 

862 and the derivation of its phase diagram triggered an explosion of research on exactly 

863 solvable models in NESM. The fundamental observation [118) is the existence of 

864 recursion relations for the stationary probabilities between systems of different sizes. 

865 These recursions are particularly striking when a = /3 = 1 [118 ; they can be 

866 generalized to arbitrary values of a and /? |12H I120j and also to the more general 

867 case in which backward jumps are allowed (also known as a partially asymmetric 

868 exclusion process, or PASEP). The most elegant and efficient way to encode these 

869 recursions is to use the matrix Ansatz |120| . We caution readers that the matrices 

870 to be presented in this subsection have nothing to do with the transition matrices M 

871 discussed above. They do not represent the intrinsic physics of TASEP, but instead, 

872 provide a convenient framework for representing the algebra that arises from the 

873 recursion relations. In particular, we need only two matrices (D and E) and two 

874 'eigenvectors' (a bra {W\ and a ket \V), normalized by (H^|V^) — 1) here. The algebra 

875 we need is 



877 
878 



DE = D + E 

(VKIE = -{W\ . (19) 
a 

876 We emphasize that in general the operators D and E do not commute with each other, 
so that an explicit representation would be typically infinite dimensional. 

Remarkably, it was shown that the stationary P* (C) for TASEP can be written 
879 as the scalar 

1 ^ 

P*{al,...,aL) = —m\\{(^^'D^{\-ai)^)\V). (20) 

r=i 

881 where 

= (W|(D + E)^|T/) (21) 

883 is a normalization constant. Each operation of the matrices D or E is associated with 

884 a filled or empty site on the lattice. For example, the weight of the configuration shown 

885 in figurelBJc) is (iy|EDDEDEDEDE|y)/Zi°. To obtain explicit values for one 

886 method is to find an explicit representation of this algebra. Another possibility is 

887 to use systematically the algebraic rules (|19p without referring to any representation. 

888 Indeed, a product of D's and E's in an arbitrary order can be decomposed as a linear 

889 combination of monomials of the type E"D™ by using repeatedly the rule DE = D+E. 

890 Then, using (iy|E"D™ \V) = a""/?""*, we find that any matrix element of the type 

891 (I20p is a polynomial in 1/a and 1//3. In particular, we find the following general 
formula 120 for Zl- 



892 



Zl = {W\{B + Ef\V) 



E 



p{2L-l-p)\p- 



Ll(L-p)l 13-^ -a- 
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894 When a — P — 1, Zi^isa. Catalan number ^120 . Another special case is a = 1 — /3, 

895 for which the scalar representation, D = 1//3,E = l/a, suffices and P* simplifies 

896 greatly. In all cases, quantities of physical interest (current, profile, correlations) can 

897 be explicitly computed using the algebra In this sense, the TASEP is 'solved:' 

898 All equal time correlations in the stationary state are known. 

899 The matrix method may look puzzling at first sight. One informal way to motivate 

900 it is the following: the steady state weight of a configuration cannot be expressed in 

901 general as a product of single site occupancy or vacancy probabilities because in the 

902 presence of interactions there are non-zero correlations (i.e., mean- field theory is not 

903 exact). However, by writing the stationary measure as a product of matrices, a sort 

904 of factorization still holds and at the same time correlations are taken into account by 

905 the fact that the matrices do not commute. 

906 3.5. Phase diagram of the open TASEP 

907 Thanks to the matrix representation method, exact expressions for the phase diagram 

908 of the TASEP, as well as stationary equal-time correlations and density profiles, can 

909 be readily obtained. For example, the mean occupation of a site i (with 1 < i < L) is 

910 given by 

(a.) ^-^{W\{T> + D (D + E)^-* \V) . (23) 

912 This expression is obtained by summing over all configurations in which site i is 

913 occupied. 

914 Similarly, the average value of the steady state current J through the "bond" 

915 connecting site i and i + 1 can be calculated as 

J{a,f3,L) = {(Ti{l - CTi+i)) 

= ^{W\ (D + Ef-' DE (D + Ef—' \V) (24) 

^ Zl-1 
Zl ' 

917 where we have used the algebra (fT9|) . We note that this value does not depend on the 

918 specific bond considered. This was expected because particles are neither created nor 

919 destroyed in the system. 

920 

921 In the limit of large system sizes, the phase diagram (figure [H consists of three 

922 main regions 

923 • In the Low-Density Phase, when a < min{/?, 1/2}, the bulk-density is p = a and 

924 the average current J = a{l — a) is a function only of the rate-limiting injection 

925 process. 

926 • In the High Density phase, when /3 < min{a, 1/2}, the typical density is 

927 characterized hy p = 1 — /3 and the steady-state current, J = /3(1 — /3), is a 

928 function only of the rate-limiting extraction step. 

929 • In the Maximal Current Phase, when a > 1/2 and /3 > 1/2, the bulk behavior is 

930 independent of the boundary conditions and one finds p = 1/2 and J = 1/4. In 

931 this phase, particle-particle correlations decay only algebraically away from the 
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932 boundaries, in contrast to exponentially decaying correlations in both the low and 

933 high density phases. 

934 • The low and high density phases are separated by the 'shock-line', a = P < 1/2, 

935 across which the bulk-density is discontinuous. In fact, the profile on this line is 

936 a mixed-state of shock-profiles interpolating between the lower density p — a and 

937 the higher density p ~ 1 ~ /3. 

938 Detailed properties of the phase diagram are reviewed in [Ml 11221 1123) . 




Figure 2. (a) Exact solution (derived from eqs. I)24| l and l|22| l for the steady-state 
current J of a TASEP with L = 10 sites, (b) The steady-state current J plotted 
in the L oo limit, (c) The phase diagram for the current of an infinite lattice 
{L = oo) TASEP as a function of the injection and extraction rates. 



939 The particle density profiles in the large L limit, described in each phase above, 

940 are only approximately uniform in that they are asymptotically accurate only in the 

941 bulk, but must vary near the boundaries in order to satisfy conditions determined by 

942 the steady-state particle injection and extraction |124| . It turns out that the MFT 

943 approaches, recursion and hydrodynamic equations, reproduce the exact L ^ oo phase 

944 diagram; however, the matrix product approach finds the correct density profile which 

945 is not obtained by mean-field approximations. 



946 3.6. Some extensions of the matrix Ansatz 

947 Extension to the general ASEP model: 

948 If we allow jumps in both directions and entrance/exit at both ends of the system, 

949 the matrix technique can still be applied. The algebra (|19|) must be replaced by the 

950 more general rules 



pDE gED = D + E 

(/3D-,5E)|y) =\V) 

{W\ (aE-7D) = {W\ . (25) 

951 This new algebra allows one to calculate the general phase diagram of the open ASEP 

952 using orthogonal polynomials |125[ 11231 11261 1127) . The phase diagram of the ASEP 

953 turns out to be identical to that of the TASEP in a two-dimensional slice across 

954 effective parameters that are functions of the intrinsic rates a, /3, 7, 5, p, g. 
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955 Multispecies exclusion processes: 

956 The matrix method can be used to represent the stationary measure of periodic 

957 systems with defects or with different classes of particles |1061 11071 11081 11091 IllOi 

958 111111112] , In particular, this allows us to investigate shock profiles that appear in the 

959 ASEP and to prove that these shocks, predicted by Burgers' equation, do exist at the 

960 microscopic level and are not artifacts of the hydrodynamic limit jl06[ I128| I129| 1130) . 

961 

962 Macroscopic density profiles: 

963 Knowing exactly the microscopic invariant measure allows us to rigorously coarse- 

964 grain the problem and to determine the probability of occurrence of a density profile 

965 that differs from the average profile. The calculation of such 'free energy functionals' 

966 is an important step in understanding how non-equilibrium macroscopic behavior 

967 emerges from the microscopic model. For a review of these very important issues 

968 and for relevant references, we refer the reader to [77l 1131) . 

969 Other models: 

970 We emphasize that the matrix product representation method has proved to be 

971 very fruitful for solving many one-dimensional systems; a very thorough review of this 

972 method can be found in [123) . 

973 3.7. Time- dependent properties: the Bethe Ansatz 

974 In order to investigate the behavior of the system away from stationarity, the spectrum 

975 of the Markov matrix is needed. For an arbitrary stochastic system, the evolution 

976 operator cannot be diagonalized analytically for any system sizes. However, the ASEP 

977 belongs to a very special class of models: it is integrable and it can be solved using 

978 the Bethe Ansatz as first noticed by D. Dhar in 1987 [132 . Indeed, the Markov 

979 matrix that encodes the stochastic dynamics of the ASEP can be rewritten in terms 

980 of Pauli matrices; in the absence of a driving field, the symmetric exclusion process 

981 can be mapped exactly into the Heisenberg spin chain. The asymmetry due to a 

982 non-zero external driving field breaks the left /right symmetry and the ASEP becomes 

983 equivalent to a non-Hermitian spin chain of the XXZ type with boundary terms that 

984 preserve the integrable character of the model. The ASEP can also be mapped into a 

985 six vertex model [551 11331 1134) . These mappings suggest the use of the Bethe Ansatz 

986 to derive spectral information about the evolution operator, such as the spectral gap 

987 [131 [1351 [nZl dSl] and large deviation functions [MOl [141] . 

988 Here, we will apply the Bethe Ansatz to the ASEP on a ring. A configuration can 

989 also be characterized by the positions of the N particles on the ring, (xi, X2, . . . ,xn) 

990 with 1<xi<X2<...<xn<L. With this representation, the eigenvalue equation 

991 (fT5]l becomes 

EiI){xi,...,xn) = 

992 X]iP[V'(a;i,---,a;i_i, - 1, a;j+i,...,XAr) -■(/'(a;i,...,x„)] -I- (26) 

9 [V'Ca:^!, • ■ • Xj + l, Xj+i, . . . ,xn) - i^ixi, . . . ,Xn)] , 

993 where the sum runs over the indices i such that Xi^i < Xi — 1 and over the indices 

994 j such that Xj + 1 < Xj+i ; these conditions ensure that the corresponding jumps are 

995 allowed. 
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996 We observe that equation (j26p is akin to a discrete Laplacian on a iV-dimensional 

997 lattice: the major difference is that the terms corresponding to forbidden jumps are 

998 absent. Nevertheless, this suggests that a trial solution, or Ansatz^ in the form of 

999 plane waves may be useful. This is precisely the idea underlying the Bethe Ansatz, 

1000 originally developed to study the Heisenberg spin chain model of quantum magnetism 

1001 jl42j . Following Bethe's method, we will also refer to ip as a. 'wavefunction,' but 

1002 the reader should not confuse our problem with one in quantum mechanics. In our 

1003 opinion, the ASEP is the one of the simplest system with which one can learn the 

1004 Bcthc Ansatz. The next subsection is devoted to such a tutorial. 

1005 3.8. Bethe Ansatz for ASEP: a crash-course 

1006 Our aim is to solve the linear eigenvalue problem (I26p which corresponds to the 

1007 relaxation modes of the ASEP with TV particles on a ring of L sites. We will study 

1008 some special cases with small N to illustrate the general structure of the solution. 

1009 

1010 The single particle case: For = 1, equation ((26)) reads 

1011 Eiljix) = p7p{x - I) + qijj{x + I) - (j) + q)'il){x) , (27) 

1012 with \ < X < L and where periodicity is assumed 

1013 ip^x + L) = ipix) . (28) 

1014 Equation (1271) is simply a linear recursion of order 2 that is solved as 

^{x) = Az^ + Bzl , (29) 

1015 where r = z± are the two roots of the characteristic equation 

qr^ -{E + p + q)r+p = 0. (30) 

1016 The periodicity condition imposes that at least one of the two characteristic values is 

1017 a i-th root of unity (Note that because z^z_ — p/q, only one of them can be a root 

1018 of unity when p ^ q). The general solution is given by 



iP(x) = Az^ with 2-^ = 1, (31) 

1019 i.e., simple plane waves with 'momenta' being integer multiples of lix jh and associated 

1020 with eigenvalue 

1021 = - + 92 - (p + 9) . (32) 

z 

1022 

1023 The two-particle case: The case N = 2 where two particles are present is more 

1024 interesting because when the particles are located on adjacent sites the exclusion effect 

1025 plays a role. Indeed, the general eigenvalue equation (j26p can be split into two different 

1026 cases: 

1027 • In the generic case, Xi and X2 are separated by at least one empty site: 

Elp{xi,X2) =p[i'{xi - l,X2) + 1p{xi,X2 - 1)] 

1028 (33) 

+q[ip{xi + 1,2:2) + ip{xi,X2 + 1)] - 2{p + q)ip{xi,X2) . 
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1029 • In the (special) adjacency case, X2 — Xi + 1, some jumps are forbidden and the 
eigenvalue equation reduces to: 



1030 



Eij^lxi.Xi + 1) = pi^ixi — 1, xi + 1) + qil){xi,xi + 2) — (p + q)il){xi,xi + 1). (34) 

This equation differs from the generic equation (j33p in which we substitute 
X2 = xi -\- 1. An equivalent way to take into account the adjacency case is 
to impose that the generic equation (j33]) is valid for all values of xi and X2 and 
add to it the following cancellation boundary condition: 

pipixijXi) + qtl^ixi + I, Xi + 1) — {p + q)'ip{xi,xi + 1) = . (35) 

We now examine how these equations can be solved. In the generic case particles 
behave totally independently (i.e., they do not interact). The solution of the generic 
equation ([55]) can therefore be written as a product of plane waves tp{xi,X2) = 
ydzJ'^Zj^, with the eigenvalue 

E^p(— + —]+q{zi+Z2)~2{p + q). (36) 

V^i Z2J 

However, the simple product solution cannot be the full answer: the cancellation 
condition for the adjacency case ([55]) also has to be satisfied. The first crucial 
observation, following H. Bethe |142) . is that the eigenvalue E, given in ((36)) is invariant 
1045 by the permutation zi ^ Z2. In other words, there are two plane waves Az^^Z2^ 
and _Bz2^zJ^^ with the same eigenvalue E which has a two-fold degeneracy. The full 
eigenfunction corresponding to E can thus be written as 

V^(xi,X2)=Ai2Zr^2'+^21^r^l% (37) 

where the amplitudes A12 and A21 are yet arbitrary. The second key step is 
to understand that these amplitudes can now be chosen to fulfill the adjacency 
cancellation condition: substitution of p7l) into equation psp . we obtain 

A21 ^ qziZ2 - {p + q)z2+p ^gg^ 
A12 qziZ2 - (p + q)zi + p 

The eigenfunction (jST]) is therefore determined, but for an overall multiplicative 
constant. We now implement the periodicity condition that takes into account the 
fact that the system is defined on a ring. This constraint can be written as follows for 

1 < Xi < X2 < L 

ip{xi,X2) ^ ip{x2,xi + L) . (39) 

This relation plays the role of a quantification condition for the scalars zi and Z2, 
which we will call Bethe roots. If we now impose the condition that the expression (|37| 
satisfies equation p9p for all generic values of the positions xi and X2, new relations 
between the amplitudes arise: 

4^ = 4 = 4 . (40) 

1063 Comparing equations ([551) and ([^(H) leads to a set of algebraic equations obeyed by the 

1064 Bethe roots zi and Z2- 

^ _ qziZ2 - {p^q)zi -\-p 



1032 
1033 
1034 
1035 

1036 

1037 
1038 
1039 
1040 



1042 
1043 
1044 



1046 
1047 

1048 

1049 
1050 

1051 



1053 
1054 
1055 
1056 

1057 

1058 
1059 
1060 
1061 



qziZ2- {p^q)z2^p 
qziZ2- {p + q)z2+p 
qziZ2 - (p + q)zi +p 



(41) 
(42) 
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1065 These equations are known as the Bethe Ansatz Equations. Finding the spectrum 

1066 of the matrix M for two particles on a ring of size L is reduced to solving these two 

1067 coupled polynomial equations of degree of order L with unknowns zi and Z2 ■ Surely, 

1068 this still remains a very challenging task but the Bethe equations are explicit and very 

1069 symmetric. Besides, we emphasize that the size of the matrix M (and the degree of 

1070 its characteristic polynomial) is of order . 

1071 

1072 The three-particle case: We are now ready to consider the case iV = 3. For a 

1073 system containing three particles, located at xi < X2 < x^, the generic equation, valid 

1074 when the particles are well separated, can readily be written using equation (|26p . But 

1075 now, the special adjacency cases are more complicated: 

1076 (i) Two particles are next to each other and the third one is far apart; such a setting 

1077 is called a 2-body collision and the boundary condition that results is identical to the 

1078 one obtained for the case N = 2. There are now two equations that correspond to the 

1079 cases xi = a; < a;2 = a; + 1 <C X3 and xi <ti X2 — x < x^ — x + 1: 

ptpix, X, X3) + q^jj{x + 1, a; + 1, 2:3) — {p + q)i^{x, x + 1, X3) — (43) 



pip{xi,x, x) + q^/j{xi,x + l,x + I) — {p + q)^{xi,x, a; + 1) 0. (44) 

1080 We emphasize again that these equations are identical to equation psp because the 

1081 third particle, located far apart, is simply a spectator {x^ is a spectator in the first 

1082 equation; xi in the second one). 

1083 (ii) There can be 3-body collisions, in which the three particles are adjacent, with 

1084 xi = x^X2 — X + 1, x:i — X + 2. The resulting boundary condition is then given by 

p [tp{x, x,x + 2)+ ij{x, x + l,x + l)]+q [tjjix + 1, a; + 1, x + 2) + ■0(a;, a; + 2, x + 2)] 

-2{p + q)^P{x,x + l,x + 2) = 0. (45) 

1085 The fundamental remark is that 3-body collisions do not lead to an independent 

1086 constraint. Indeed, equation (1451) is simply a linear combination of the constraints (|43p 

1087 and p4)) imposed by the 2-body collisions. To be precise: equation (|45|) is the sum 

1088 of equation (|43p. with the substitutions x ^ x and X3 — > a; + 2, and of equation ()44|) 

1089 with xi ^ X and a: — ^ a; + 1. Therefore, it is sufficient to fulfill the 2-body constraints 

1090 because then the 3-body conditions are automatically satisfied. The fact that 3-body 

1091 collisions decompose or 'factorise' into 2-body collisions is the crucial property that 

1092 lies at the very heart of the Bethe Ansatz. If it were not true, the ASEP would not 

1093 be exactly solvable or 'integrable'. 

1094 For = 3, the plane wave ip{xi,X2, x^) = Az^^ Z2^ z^^ is a solution of the generic 

1095 equation with the eigenvalue 

E=p(— + ^—^]+q{zi+Z2 + Z3)-3{p + q). (46) 

\Zi Z2 + Z3 / 

1097 However, such a single plane wave does not satisfy the boundary conditions (j43p 

1098 and (|44p . Again, we note that the eigenvalue E is invariant under the permutations 

1099 of zi,Z2 and Z3. There are 6 such permutations, that belong to S3, the permutation 

1100 group of 3 objects. The Bethe wave-function is therefore written as a sum of the 6 

1101 plane waves, corresponding to the same eigenvalue E, with unknown amplitudes: 
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1102 
1103 
1104 
1105 
1106 
1107 
1108 

1109 

1110 
1111 
1112 
1113 
1114 
1115 
1116 



1117 
1118 
1119 



1123 



1128 
1129 
1130 
1131 



/ ™ ^ \ A ^1 ^2 X-'i t A ^1 ^2 ^3 \ A Xi Xo Xr 

^p(xi,X2,X3) = Ai23 z^^z^'z^^ + Ai32 z/zg^Za'' + A213 Z2 -^1 -^a 

I \ ^Xi X2 X3 , A ^Xi X2 X3 , A „Xi X2 X3 

+ A231 Z2 Z3 Zi + ^312 2^3 Zi Z2 + A32I Z3 Z2 Zi 

seS3 

The 6 amplitudes As are uniquely and unambiguously determined (up to an overall 
multiplicative constant) by the 2-body collision constraints. It is therefore absolutely 
crucial that 3-body collisions do not bring additional independent constraints that 
the Bcthe wave function could not satisfy. We encourage the reader to perform the 
calculations (which are very similar to the N = 2 case) of the amplitude ratios. 

Finally, the Bethe roots zi, Z2 and Z3 are quantized through the periodicity 
condition 

lp{xi,X2,X3) ^ lp{x2,X3,Xi + L) , (48) 

for 1 < xi < X2 < X3 < L. This condition leads to the Bethe Ansatz equations (the 
equations for general N are given below). 

The general case: Finally, we briefly discuss the general case iV > 3. Here one 
can have /c-body collisions with k — 2,3,...N. However, all multi-body collisions 
'factorize' into 2-body collisions and ASEP can be diagonalized using the Bethe Wave 
Function 

■lp{xi,X2, . . . , Xn) ^ ^ As Zslifs{2) ' ' ' ^s(Ar) ' (^9) 
sGSjv 

where Sn is the permutation group of N objects. The iV! amplitudes As are fixed 
(up to an overall multiplicative constant) by the 2-body collisions constraints. The 
corresponding eigenvalue is given by 

N ^ N 

E = pY,-+qY,z,-N{p + q). (50) 
1=1 * 1=1 
The periodicity condition 

V'(xi,a;2, . . ■,xn) = ip{x2,X3, . . . ,xn,xi + L) , (51) 



where 1 < xi < X2 < ■ ■ ■ < Xn < L, leads to a set of algebraic equations satisfied by 
1124 the Bethe roots zi, Z2, . . . , z^r. The Bethe Ansatz equations are given by 

i-},qz^Zj-{p + q)zj+p' 

1126 for i = 1, . . . A^. The Bethe Ansatz thus provides us with a set of N coupled algebraic 

1127 equations of degree of order L (Recall that the size of the matrix M is of order 2^, 
when N ~ L/2). Although the degree of the polynomials involved are extremely high, 
a large variety of methods have been developed to analyze them j88j I143[ 1144) . 

We remark that ior p = q = 1 the Bethe equations are the same as the ones derived 
by H. Bethe in 1931. Indeed, the symmetric exclusion process is mathematically 
1132 equivalent to the isotropic Heisenberg spin chain (although the two systems describe 
totally different physical settings). 



1133 
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1143 



1145 
1146 
1147 
1148 
1149 
1150 



1152 
1153 
1154 



1160 

1161 



3.9. Some applications of the Bethe Ansatz 

The Bethe Ansatz allows us to derive information about the spectrum of the Markov 
matrix that governs the evolution of ASEP. Below, we review some applications and 
refer to the original works. 

Structure of the Bethe wave function: 

In the totally asymmetric case (TASEP), the Bethe equations simplify and it is 
possible to perform analytical calculations even for finite values of L and TV. Besides, 
the TASEP Bethe wave function takes the form of a determinant: 

V'(ei,...,eJv) = det(R), (53) 

where R is a x matrix with elements 

^^^ = 7T^^ forl<z,j<7V, (54) 

(zi, . . . , zn) being the Bethe roots. By expanding the determinant, one recovers the 
generic form ((49|) for the Bethe wave function. The formulas (|53|) and (|54[) can be 
verified directly by proving that the eigenvalue equation (1261) and all the boundary 
conditions are satisfied [11211145] . As we will see in the next subsection, determinantal 
representations of the eigenvectors and of the exact probability distribution at finite 
time play a very important role in the study of the TASEP (HSl [Wl [TIHlfng] . 



1151 Spectral gap and dynamical exponent: 

The Bethe Ansatz allows us to determine the spectral gap which amounts to 
calculating the eigenvalue Ei with the largest real part. For a density p — N /L, one 
obtains for the TASEP 



-6.509189337...^ 2i7r(2p- 1) 



E,=-2^^^^)- -^j^—l±^Ji L. (55) 

^ ^ ^ ^ ^ ■' 

relaxation oscillations 

1156 The first excited state consists of a pair of conjugate complex numbers when p is 

1157 different from 1/2. The real part of Ei describes the relaxation towards the stationary 

1158 state: we find that the largest relaxation time scales as T ~ with the dynamical 

1159 exponent z = 3/2 [1501 11321 1135| I136[ I138[ . This value agrees with the dynamical 
exponent of the one-dimensional Kardar-Parisi-Zhang equation that belongs to the 
same universality class as ASEP (see the review of Halpin-Healy and Zhang 1995 

1162 |l5l)). The imaginary part of Ei represents the relaxation oscillations and scales as 

1163 L~^] these oscillations correspond to a kinematic wave that propagates with the group 

1164 velocity 2p—l\ such traveling waves can be probed by dynamical correlations 152 . 99] . 

1165 The same procedure also allows us to classify the higher excitations in the spectrum 

1166 [153J • For the partially asymmetric case {j>,q > 0,p ^ g), the Bethe equations do not 

1167 decouple and analytical results are much harder to obtain. A systematic procedure 

1168 for calculating the finite size corrections of the upper region of the ASEP spectrum 

1169 was developed by Doochul Kim [136| 1137] . 

1170 Large deviations of the current: 

1171 The exclusion constraint modifies the transport properties through the ASEP 

1172 system. For instance, if one considers the ASEP on a ring and tags an individual 

1173 particle (without changing its dynamics), the particle displays diffusive behavior in 

1174 the long time limit, but with a tracer diffusion coefficient 2? that depends on the 
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size of the system. It is different from the free, non-interacting particle diffusion 
constant: I?asep 7^ 'Dircc- In a ring of size L, I?asep scales as in presence of 

an asymmetric driving. In the case of symmetric exclusion, one has Psep ~ L~^. 
These scaling behaviors indicate that in the limit of an infinite system the diffusion 
constant vanishes and fluctuations become anomalous jl54) . In the classic result 
[91] [92I l97] of the SEP on an infinite line, root-mean square displacement of a tagged 
particle scales as t^^^. In the ASEP, subtleties associated with the initial conditions 
arise. For a fixed initial condition, the position of the tagged particle spreads as t^^^ . 
However, if an average is carried out by choosing random initial conditions (with 
respect to the stationary measure), then a normal i^/'^ diffusion is recovered. The 
latter is a trivial effect due to the local density fluctuations that result from varying 
the initial condition; these fluctuations completely overwhelm and mask the intrinsic 
dynamical fluctuations. Another feature that one expects is a non-Gaussian behavior, 
i.e., cumulants beyond the second are present. 

An alternative observable that carries equivalent information (and is amenable to 
analytical studies) is the total current transported through the system. Consider the 
statistics of Yi, the total distance covered by all the particles between the time and 
t. It can be shown (see e.g., jll2j and references therein) that in the long time limit 
the cumulant generating function of the random variable Yj behaves as 

(e^^')~e^(^)*. (56) 

This equation implies that all the cumulants of Yf grow linearly with time and can be 
determined by taking successive derivatives of the function E{^) at /i = 0. Another 
way to characterize the statistics of Yt is to consider the large- deviation Junction of 
the current, defined as follows 

Prob (^=j^^ e-*^(j) . (57) 

From equations (I56l57p . we see that G{j) and i?(/i) are the Legendre transforms of 
each other. Large deviation functions play an increasingly important role in non- 
equilibrium statistical physics 78 , especially through application of the Fluctuation 
Theorem 76J. Thus, determining exact expressions for these large deviation functions 
for interacting particle processes, either analytically or numerically, is a major 
challenge in the field [77|. Moreover, higher cumulants of the current and large 
deviations are also of experimental interest in relation to counting statistics in quantum 
systems !158j. 

The crucial step in the calculation of the cumulants is to identify the generating 
function E{fi) as the dominant eigenvalue of a matrix M(/^), obtained by the following 
deformation of the original Markov matrix M: 

M(/x) = Mo + e^Mi + e"^M_i, (58) 

where Mq is the matrix of the diagonal of M, and Mi (M_i) is a matrix containing the 
transitions rates of particle hopping in the positive (negative) direction. Hence, the 
statistics of the transported mass has been transformed into an eigenvalue problem. 
The deformed matrix M(/i) can be diagonalized by the Bethe Ansatz by solving the 
following Bethe Ansatz equations 

^J-^ qe >^ZiZj - (p + q)Zj + pe'^ 
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1225 
1226 



1228 
1229 

1230 

1231 
1232 



1246 



The corresponding eigenvalues are given by 

N N 



E{ii;zi,Z2...ZN)=pe^'^ — + qe^'^z,-N{p + q). (60) 

i=l * i=l 

1220 For = we know that the largest eigenvalue is 0. For /i 7^ the cumulant generating 

1221 function corresponds to the continuation of this largest eigenvalue (the existence of 

1222 this continuation, at least for small values of /i is guaranteed by the Perron-Frobenius 

1223 theorem) . 

1224 We remark that equations (j59p and (j60p are invariant under the transformation 

1 

z 

IJ ^ fiQ — with /iQ = log - . (61) 

P 

This symmetry implies that the spectra of M(/i) and of M(/io — m) are identical. This 
functional identity is in particular satisfied by the largest eigenvalue of M and we have 

Eifi) = E{^io - m) ■ (62) 

Using the fact that the LDF G{j) is the Legendre transform of E{fj.), we obtain the 
canonical form of the Fluctuation Theorem (or Gallavotti-Cohen relation) 

G{j) = G{-j) ~ fioj ■ (63) 

We observe that in the present context this relation manifests itself as a symmetry of 
the Bethe equations. However, it is satisfied by a large class of systems far from 

1233 equilibrium [7TJ [751 [73]. The validity of the Fluctuation Theorem does not rely 

1234 on the integrability of the system and it can be proved for Markovian or Langevin 

1235 dynamical systems using time- reversal symmetry [75[ 176) as well as measure-theoretic 

1236 considerations jl59( 1160] . 

1237 The complete calculation of the current statistics in the ASEP is a difficult 

1238 problem that required more than a decade of effort. The breakthrough was the solution 

1239 of the TASEP, by B. Dcrrida and J. L. Lcbowitz in 1998 [140 . These authors obtained 

1240 the following parametric representation of the function E{fi) in terms of an auxiliary 

1241 parameter B: 

fe=i ^ ^ 
1243 where B{fi) is implicitly defined through 

^--E(fc;^)|^- (65) 

k=l ^ ^ 



1245 These expressions allow us to calculate the cumulants of Yt, for example the mean- 
current J and a 'diffusion constant' D: 



J ^ lini (^*) _ d£;(/i) _ N{L - N) 



t-i-oo t d/x 
D = lim 



(66) 



t d/i2 



M=o (L-l)!2 (2iV-l)! (2L-2iV-l) 



NESM: A paradigm and applications 



32 



1258 
1259 
1260 
1261 
1262 
1263 
1264 
1265 



When L ^ oo, p — L/N, and \ j ~ Lp{l — p)\ <^ L, the large deviation function G{j) 
can be written in the following scaling form: 

GU) - xf^Hfl^ML^) (67) 



ttN^ \ p(l - p) 



with 



2 

i/(y)^-^|y|3/2 for y^-oo. (69) 

1252 The shape of the large deviation function is skewed: it decays with a 5/2 power as 

1253 y — >■ +00 and with a 3/2 power as j/ — ?■ — oo. For the general case 9 7^ on a ring, 

1254 the solution was found by rewriting the Bethe Ansatz as a functional equation and 

1255 restating it as a purely algebraic problem |16HI162] . For example, this method allows 

1256 us to calculate the first two cumulants, J and D: 



J=(^^^(^^fc^Lp(l-p)for L 
p L — \ p 



D = 



fc>0 ^ ^ 



(70) 



where are the binomial coefficients. Note that the formula for D was previously 
derived using the matrix product representation discussed above [1631 . Higher 
cumulants can also be found and their scaling behavior investigated. For ASEP 
on a ring, the problem was completely solved recently by S. Prolhac [TMj . who 
found a parametric representation analogous to equations (1641) but in which the 
binomial coefficients are replaced by combinatorial expressions enumerating some tree 
structures. A particularly interesting case |165) is a weakly driven limit defined by 
q/p— 1 — j^,L—i'Oo. In this case, we need to rescale the fugacity parameter as /i/L 

1266 and the following asymptotic formula for the cumulant generating function can be 

1267 derived 



^ IZ' ^ - ZJ " L 2l^ + - + ' 

00 

and where -Bj's are Bernoulli Numbers. We observe that the leading order (in 1/L) 
is quadratic in fi and describes Gaussian fiuctuations. It is only in the subleading 
correction (in 1/i^) that the non-Gaussian character arises. This formula was also 
obtained for the symmetric exclusion case = in |:166j . We observe that the series 
that defines the function (p{z) has a finite radius of convergence and that (p{z) has a 
singularity for z = — vr^. Thus, non-analyticities appear in E(p, v) as soon as 

27r 



- P) 



1268 By Legendre transform, non-analyticities also occur in the large deviation function 

1269 G{j) defined in ([57)1 . At half-filling, the singularity appears at Vc = 47r as can be seen 
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in figure |31 For v < the leading behavior of G{j) is quadratic (corresponding to 
Gaussian fluctuations) and is given by 



ALp{l - p) 



(73) 



1273 For V > Vc, the series expansions (|71I72I) break down and the LDF G{j) becomes 

1274 non-quadratic even at leading order. This phase transition was predicted by T. 

1275 Bodineau and B. Derrida using macroscopic fluctuation theory |167[ 1168] I169j . These 

1276 authors studied the optimal density profile that corresponds to the total current j 

1277 over a large time t. They found that this optimal profile is flat for j < jc but it 

1278 becomes linearly unstable for j > jc- In fact, when j > jc the optimal profile is 

1279 time-dependent. The critical value of the total current for which this phase-transition 

1280 occurs is jc = p(\ — p) ^ — and therefore one must have > i^c for this transition 

1281 to occur. One can observe in figure [3] that for v ~> the large deviation function 

1282 G'(j) becomes non-quadratic and develops a kink at a special value of the total current 

1283 J . 




7>fvp(l-p)) 



7/(vp(l-p)) 



7Xvp(l-p)) 



Figure 3. Behaviour of the large deviation function G as a function of the current 
i I (iyp{l — p)) for different values of u. The gray dots correspond to L = 50, N = 25 
and the black dots correspond to L = 100, N = 50. The curves are formally 
obtained by numerically solving the Bethe Ansatz equations I I59I1 . and Legendre 
transforming E(fi). The thin blue curve represents the leading Gaussian behavior 
indicated by equation l|73| l. 



1284 Bethe Ansatz for other systems out of equilibrium: 

1285 We have discussed in this review the Bethe Ansatz for a system on a periodic ring 

1286 where the total number of particles is a conserved quantity. It is possible to extend 

1287 the Bethe Ansatz for the finite ASEP with open boundaries in which the total number 

1288 of particles is not constant. The resulting Bethe equations have a more complicated 

1289 structure; they have been thoroughly analyzed in a series of papers by J. de Gier and 

1290 F. Essler |170[ 11531 1171j who calculated the spectral gaps in the different phases of 

1291 the open ASEP. In particular, they discovered sub-phases in the phase diagram that 

1292 could not be predicted from the steady state measure alone. We note that the Bethe 

1293 Ansatz can be applied to variants of the ASEP, such as models with impurities |141| . 

1294 multispecies exclusion processes |1721 11731 1174| as the zero-range process [175] , the 

1295 raise and peel model |125j . vicious walkers [176', or the Bernoulli matching model of 

1296 sequence alignment |177i 1178] . 
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1297 3.10. ASEP on an infinite lattice: Bethe Ansatz and random matrices 

1298 In this last subsection, we briefly review some important analytical results that have 

1299 been obtained for exclusion processes on an infinite line, especially in connection with 

1300 the Bethe Ansatz method discussed above, this approach will allow us to derive 

1301 some insight into the dyanmics of the systems. More detailed results and historical 

1302 discussions can be found in the literature cited in the text. 

1303 We first consider the case of the ASEP on the infinite lattice Z but with only a 

1304 finite number TV of particles. Because the particles cannot overtake one another, the 

1305 ordering of the particles is conserved by the dynamics, and we can label the particles 

1306 from right to left. Suppose that at t = 0, the particles are located at positions 

1307 yi > . . . > j/jv and at some later time they are located at xi > . . . > xjv- 

1308 The transition probability, P(a;i, . . . , xjv, . • . , J/at, 0), for reaching the final 

1309 configuration . . . , xm at t starting from . . . , y^v at i = has the following exact 

1310 expression: 



1311 The crucial observation is that there exists a closed formula for the amplitude A^, 

1312 given by 



1328 
1329 
1330 
1331 
1332 
1333 
1334 



As{{z}) ^ sgn s — , (75) 

1313 where we use the convention p + g = 1 . The expressions (|74p and (j75p are reminiscent 
of the formulae given by the Bethe Ansatz. These results were initially derived for 
the TASEP on Z by G. Schiitz in 1997, who constructed them from the Bethe Ansatz 
jl47| and then proved rigorously that equation ([71]) solves exactly the time-dependent 
Markov equation with the correct initial condition (generalized to the periodic case 
by V. Priezzhev in |148[ 1149] ). The general result for the ASEP was found by Tracy 
and Widom in 2008 [179] and has motivated many studies in the last few years. 
Equation (|74p is a seminal result: it is an exact expression for the Green function of 
the ASEP from which, in principle, individual particle distributions and correlations 
functions can be extracted (this can be an extremely difficult task in practice) . Finally, 

1323 we emphasize that the 'Bethe roots' zi are not quantified in the infinite system: each 

1324 of them takes a continuous set of values along the circular contour of radius i?, so 

1325 small that the poles of As({z}) lie outside C^. 

1326 To give one specific example, consider the totally asymmetric exclusion with unit 

1327 hopping rate p = 1 (and q = 0) and with a step initial condition, where all sites 
right of the origin {i > 0) are empty and all to the left {i < 0) are occupied. If we 
are interested in the behaviour of only the right most N particles after time t, then 
we can replace the semi-infinite string by a finite segment of N particles. The point 
is, in a TASEP, none of the particles hops to the left and so, particles to the left of 
our A'^-particle string cannot affect their behaviour. Thus, it is sufficient to consider 
P{xi, . . . , xn, t\yi = 0, . . . ,yN = 1 — N,0). Now, suppose we ask a more restricted 
question: What is P{M, N,t), the probability that the A^-th particle (initially located 



1314 
1315 
1316 
1317 
1318 
1319 
1320 
1321 
1322 



NESM: A paradigm and applications 



35 



1337 
1338 



1344 



1349 
1350 
1351 
1352 
1353 



1357 
1358 



at z = 1 — A^, has performed at least M hops by time t? In terms of P, we write 

P{M,N,t)^ J2 Pixu...,XN,t\yi=0,...,y„ = l-N,0) (76) 



a;i>...>a:iv>A/-A' 



We next exploit a determinantal representation, analogous to (l53l) and ((54|) for and 
express P as: 



1339 P{xi, . . . ,XN,t\yi, . . . ,yN,0) = det{M) , (77) 

1340 where B is a x matrix with elements 

B,j=<L -^2^'-!^^(l-z)J-*e(^-i)* for l<i,j<iV, (78) 
Jcr 27r«z 

1342 After some calculation |180) . these equations allow us to express P{M, N, t) compactly: 

N 

P{M,N,t)^ / d^x J] ix,-XjfY\xf^-^e-^\ (79) 

Zaln Jio.t]^ ,^A^^r 



l<j<j<Af j=l 



where Zm,n is a normalization factor. The integral on the right hand side has a direct 

1345 interpretation in random matrix theory: it is the probability that the largest eigenvalue 

1346 of the random matrix AA^ is less than or equal to t, with A being a. N x M matrix of 

1347 complex random variables with zero mean and variance 1/2 (unitary Laguerre random 

1348 matrix ensemble). 
Finally, we note that P{N, N, t) can also be interpreted as the probability that 

the integrated current Qt through the bond (0,1) is at least equal to A^ {Qt is the 
number of particles that have crossed the bond (0,1) between time and t). From the 
classical Tracy- Widom laws for the distribution of the largest eigenvalue of a random 
matrix [181], one can write 



t ^1/3 

where the distribution of the random variable x is given by 

Prob(x<.s) = l-^^2(-s). (81) 

The Tracy- Widom function F2{s) is the cumulative distribution of the maximum 
eigenvalue Xmax in the Gaussian Unitary Ensemble (self-adjoint Hermitian matrices) 
i.e. 

where A' ^ 1 represents here the size of the random matrix. Exact expressions 
for F2{s) can be found, for example, in |182) . This crucial relation between random 
matrix theory and the asymmetric exclusion process, first established by K. Johansson 
in 2000 |183) . has stimulated many works in statistical mechanics and probability 
theory (see e.g., [1841 11851 1186| ). The mathematical study of the infinite system has 
grown into a subfield per se that displays fascinating connections with random matrix 
theory, combinatorics and representation theory. We note that K. Johansson did 
not use the Bethe Ansatz in his original work. He studied a discrete-time version of 
the TASEP in which the trajectories of the particles were encoded in a waiting-time 
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1370 matrix, which specifies how long a given particle stays at each given site. This integer- 

1371 valued matrix can be mapped via the Robinson-Schensted-Knuth correspondence into 

1372 a Young Tableau. The value of the total current through the bond (0,1) is linearly 

1373 related to the length of the longest line of this Young Tableau. The statistics of the 

1374 length of the longest line can be found by using asymptotic analysis techniques d 

1375 la Tracy- Widom. In fact, a closely related question, the classical Ulam problem of 

1376 the longest increasing subsequence in a random permutation, was solved a few years 

1377 earlier by J. Baik, K. Deift and K. Johansson, using related methods [187j . 

1378 Very recently, in a series of articles [173 [Ml [US [IHOl [M], C. A. Tracy 

1379 and H. Widom have generalized Johansson's results to the partially asymmetric 

1380 exclusion process by deriving some integral formulas for the probability distribution 

1381 of an individual particle starting from the step initial condition. These expressions 

1382 can be rewritten as a single integral involving a Fredholm determinant that is 

1383 amenable to asymptotic analysis. In particular, a limit theorem is proved for 

1384 the total current distribution. Going from TASEP to ASEP is a crucial progress 

1385 because the weakly asymmetric process leads to a well-defined continuous limit: 

1386 the Kardar-Parisi-Zhang (KPZ) equation, a universal model for surface growth. 

1387 Indeed, a very important outgrowth of these studies is the exact solution of the 

1388 KPZ equation in one dimension. The distribution of the height of the surface at 

1389 any finite time is now known exactly (starting from some special initial conditions), 

1390 solving a problem that remained open for 25 years: a description of the historical 

1391 development of these results and the contributions of the various groups can be found 

1392 in [Ml [Hg [I93l [Ml [Ml [Ml [Ml IMl IM]- One important feature to keep in 

1393 mind is that the results (and therefore the universality class) depend strongly on the 

1394 chosen initial conditions. Lately, n-point correlation functions of the height in KPZ 

1395 have also been exactly derived by H. Spohn and S. Prolhac |2001l201j . For an overview 

1396 of these fascinating problems, we recommend the article by T. Kriecherbauer and J. 

1397 Krug [Hg, and the reviews by D. Aldous and P. Diaconis [202] and I. Corwin [203] . 

1398 3.11. Hydrodynamic mean field approach 

1399 Though elegant and rigorous, the methods presented above cannot be applied to 

1400 systems much more complex than the TASEP. Typically, further progress relies 

1401 on a very successful approach, based on the mean field approximation and a 

1402 continuum limit, leading to PDE's for various densities in the system. Known as the 

1403 hydrodynamic limit, such equations can be 'derived' from the master equation (|14j) . 

1404 The strategy starts with the exact evolution equation for pi{t) = {o'i)t = X^c '^iP {C,t). 

1405 Differentiating this Pi{t) with respect to t and using equation (|14p. we see that new 

1406 operators appear: 

0{C') = ^ai{C)M{C,C'). (83) 
c 

1408 In the case of ASEP, M{C,C') is sufficiently simple that the sum over C can be 

1409 easily performed, leaving us with products of cr's (associated with configurations 

1410 C). Applying the mean field approximation (equation (|17|) '). the right-hand side now 

1411 consists of products of Pi{t) and Pi±i{t). 

1412 Taking the thermodynamic and continuum limit, we let i, L — >■ oo with finite x = 

1413 i/L and write Pi (t) ~ p{x,t). Next, we expand pi±i(t) ~ p{x,t)zLedxP+i£'^ /2)d'^p+. . 
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1414 where e = 1/L — 0. The result is a hydrodynamic PDE for the particle density [124) : 



1416 



1418 
1419 
1420 
1421 
1422 
1423 
1424 



1427 
1428 



1430 
1431 
1432 
1433 
1434 
1435 
1436 
1437 
1438 
1439 



dp{x,t) d r 1 d 



e 



[(q-p)p(l-p)] 



dt dx '^'^^ '^'J 2 dx 



(1 - p) ^ [(p + q)p\ + {p + l)P^ 



(84) 



where e = 1/L — !■ 0. Note that even slow variations in the hopping rates {pi — >■ p{x) 
1417 and Qi <l{x)) can be straightforwardly incorporated |204| . Analogous hydrodynamic 
equations have also been derived for ASEPs with particles that occupy more than one 
lattice site [2051 12061 1207] . In addition to spatially slowly varying hopping rates, 
this approach is useful for exploring more complex models, e.g., a TASEP with 
particles that can desorb and/or attach to every site (Section 14. 2p . Of course, such 
equations are also the starting point for very successful field theoretic approaches to 
dynamic critical phenomena near-equilibrium [208[ 12091 12101 1211) as well as ASEPs 
(even with interacting particles) in higher dimensions 1^1^ WM \STE[ 

1425 Supplemented with noise terms, these become Langevin equations for the density or 

1426 mangetisation fields. Then, fiuctuations and correlations can be computed within, 
typically, a perturbative framework. 

Returning to equation (|84|) . it is especially easy to analyse in the steady-state 
1429 limit where the resulting equation is an ODE in x. But, the highest power of e 
multiplies the largest derivative, so that we are dealing with 'singularly perturbed' 
differential equations |1241 12181 1219| . The solutions for the steady-state mean field 
density p{x) consist of 'outer solutions' that hold in the bulk, matched to 'inner 
solutions' corresponding to boundary layers of thickness e at the open ends. For 
the ASEP, this approach correctly determines the phase boundaries and the self- 
consistently determined steady-state current J (which arise both in the integration 
constant of equation (j84p and in the boundary conditions). However, as anticipated, 
the density profiles are not exactly determined, even when L ^ oo. Nevertheless, for 
more realistic physical models such as the ones to be discussed in the next section, such 
a mean-field approach (or some of its variants) is often the only available analytical 
1440 method at our disposal. 



1441 4. Biological and related applications of exclusion processes 

1442 In this section, we review a number of applications of stochastic exclusion models 

1443 to problems of transport in materials science, cell biology, and biophysics. While 

1444 quantitative models of these real-world applications often require consideration 

1445 of many microscopic details (and their corresponding parameters), simple one- 

1446 dimensional lattice models nonetheless can be used to capture the dominant 

1447 mechanisms at play, providing a succint representation of the system. Moreover, these 

1448 types of models can be extended in a number of ways, and are amenable to concise 

1449 analytic solutions. As we will see, application of lattice models to complex biophysical 

1450 systems is aided by a few main extensions and modifications to the TASEP presented 

1451 above. These modified dynamics are illustrated respectively, by figures |4]|9l 

1452 (i) Longer-ranged interactions: Objects represented by particles in an exclusion 

1453 process may have molecular structure that carry longer-ranged particle-particle 

1454 interactions. For example, to model ribosomes on an mRNA or cargo-carrying 

1455 molecular motors, we should use large particles taking small discrete steps and 

1456 introduce a rule beyond on-site exclusion. We may regard these as 'extended 

1457 particles' that occlude £ > 1 lattice sites (figure |4]). 
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Figure 4. An interior section of an asymmetric exclusion process with extended 
particles of size £ = 3. The individual particles occlude three of the lattice sites 
on which the hops occur. 



(ii) Particle detachment and attachment: Particles such as molecular motors 
have finite processivity. That is, they can spontaneously detach from their lattices 
before reaching their destination. Conversely, particles in a bulk reservoir can also 
attach at random positions on the lattice. The coupling of the particle number 
to a bulk reservoir is analogous to the problem of Langmuir kinetics |219i 1220] 
on a one-dimensional lattice, except that the particles are directionally driven on 
the substrate (figure [5]). 

particle ' J vacancy ^ 

•. ^ >^ ^ ^ 



cN-C; O • O O • • 



Figure 5. A TASEP with Langmuir kinetics where particles can spontaneously 
detach and attach at every site with rates ujd and lu\, respectively. Adapted from 
[221]. 



(iii) Multiple species: Usually, biological transport involves multiple interacting 
species in confined geometries, often one-dimensional in nature. For example, 
these secondary species may represent particles that are co-transported or 
counter-transported with the primary particles, or, they may represent species 
that bind to pores and regulate primary particle transport [222] (figure [6]) . 
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Figure 6. A two-species exclusion model, adapted from |223j . where -|- and 
— particles move in opposite directions but can interconvert with rate k±. 
Attachment and detachment are also allowed with rates c± and d±, respectively. 

1470 (iv) Partial exclusion &: coupled chains: Often, the strict one-dimensional nature 

1471 of the exclusion dynamics can be relaxed to account for particles that, while 

1472 strongly repelling each other, can on occasion pass over each other. This 

1473 scenario might arise when pores are wide enough to allow particle passing. 
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Related extensions of single-chain exclusion processes are exclusion processes 
across multiple interacting chains (figure [7]). 
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Figure 7. Two coupled TASEP lattices with interchain hopping rates uii and 
W2. Adapted from |224| . 

1476 (v) Internal states and nonexponentially distributed dwell times: The 

1477 physical mechanisms of how excluding particles are inserted, extracted, and hop in 

1478 the interior of the lattice are often complex, involving many intermediate chemical 

1479 steps |225[|226j . Therefore, the distribution of times between successive hops, even 

1480 when unimpeded by exclusion interactions, is often non-exponential. Specific 

1481 hopping time distributions can be incorporated into lattice particle models by 

1482 explicitly evolving internal particle states (figure [5]) . 








i 







Figure 8. Internal states (1-5) that determine the timing between particle hops 
(adapted from 12261 V Here, a linear sequence of Poisson processes generates a 
Gamma-distributed |227| interhopping time distribution (dwell time). This model 
for the internal dynamics has been used to model mRNA translation by ribosome 
enzymes. 

1483 (vi) Variable lattice lengths: The lattices on which the exclusion processes occur 

1484 may not be fixed in some settings. A dynamically varying length can arise in 

1485 systems where growth of the lattice is coupled to the transport of particles to 

1486 the dynamically- varying end of the lattice |228l 12291 12301 1231] . While no exact 
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solutions have been found, mean-field and hydrodynamic approaches have been 
successfully applied [232] . In the continuum limit, this problem is analogous to the 
classic Stefan problem [2321 1233) , except with nonlinear particle density dynamics 
(figure ED. 



Figure 9. Schematic of an ASEP (adapted from 12320 witli detactiment and 
attachment rates fc_ and fc_|_, respectively, and a 'confining' wall at site L. The 
wall has intrinsic dynamics described by hopping rates w± . 



1491 Not surprisingly, the techniques presented in the previous section are typically not 

1492 powerful enough for obtaining exact results for these variants of the basic exclusion 

1493 processes. Of course, for sufhciently small systems (e.g., L < 10), high performance 

1494 computers can be used to diagonalize M, but that method is not viable for many 

1495 physical or biological systems. There are also some cases, including special cases of 

1496 two-species exclusion and ASEP with extended particles on a ring (to be discussed 

1497 below, within their own contexts), where exact solutions can be found for all L. Despite 

1498 the lack of progress in finding analytic results, there are two other approaches which 

1499 provide valuable insights into these more complex yet more realistic systems. One 

1500 is computer simulations, exploiting Monte Carlo techniques on our lattice models. 

1501 Once the stochastic rules of a model are specified, then, this technique corresponds 

1502 to implementing equation ([1]). In addition, molecular dynamics simulations with off- 

1503 lattice systems also proved to be very successful. The other approach is based on 

1504 designing good approximation schemes for attacking exact equations (or expressions 

1505 for specific quantities). Mean field theory (MFT), or mean field approximation, is 

1506 most often used. An example is equation (|17p for the TASEP. However, systematically 

1507 improving these approximations it is not generally straightforward. Thus, we will see 

1508 that, for TASEP with extended objects, the substitution PT)) is very poor (when 

1509 compared to simulation results). Instead, a more sophisticated scheme, implemented 

1510 by Gibbs, et.al. [13l [14], is much more successful at predicting, say, the average 

1511 current. In turn, when this scheme fails to predict other quantities, another level 

1512 of approximation |103j (in the spirit of the Bethe- Williams approximation for the 

1513 Ising model) is designed. Nevertheless, if the predictions are in good agreement with 

1514 simulations, the MFT provides some confidence that we are 'on the right track' towards 

1515 our goal: understanding these generalized ASEPs which have found wide applicability 

1516 in modeling complex biological processes. 

1517 4-^- Pore Transport 

1518 Transport of atoms and small molecules arises in both man-made structures and 

1519 cell biological systems. Materials such as zeolites form networks of one-dimensional 

1520 channels within which molecules, such as light hydrocarbons, can pass through and/or 

1521 react. A number of authors have used exclusion processes as simple descriptions of 

1522 single-file or near single- file transport. It has long been known that tracer diffusion in 
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1523 single-file channels follows subdiffusive dynamics |2341 12351 12361 12371 Wl\ . A mean- 

1524 square displacement of the form {x^) ^ ^/t \s found by considering equilibrium 

1525 fluctuations of the density around the tagged particle in an infinite system. For steady- 

1526 state particle transport across finite pores, the lattice exclusion process has proved a 

1527 more useful starting point. Given the frequent application of exclusion processes to 

1528 pore transport, two important points should be stressed. 

1529 First, the standard lattice exclusion process assumes that waiting times are 

1530 exponentially distributed and, other than lattice site exclusion, are independent. This 

1531 assumption can fail if, say, an attractive interaction between particles is comparable 

1532 to the interaction between substrate (e.g., atoms that make up the pore walls) and 

1533 the particles. Here, concerted motion can arise and has been shown to be important 

1534 in particle transport. For example, Sholl and Fichthorn have shown using molecular 

1535 dynamics (MD) simulations how concerted motions of clusters of water affects its 

1536 overall transport in molecular sieves |238j . Similarly, Sholl and Lee also showed 

1537 that concerted motions of clusters of CF4 and Xe in AIPO4-5 and AIPO4-3I zeolites, 

1538 respectively, contribute significantly to their overall mobilities |239j . Concerted motion 

1539 has also been predicted to occur in transport across carbon nanotubes |240j . These 

1540 concerted motions arise from frustration due to a mismatch between particle-particle 

1541 and particle-pore interaction ranges, allowing for a lower barrier of motion for bound 

1542 pair of particles that for an isolated, individual particle. Although concerted motion 

1543 arises from geometrically complicated arrangements of particles that form low energy 

1544 pathways in the high dimensional energy landscapes, if a small number of these 

1545 pathways support most of the probability flux, simplifying assumptions can be made. 

1546 For example, coarse-grained treatments of concerted motion were developed by Sholl 

1547 and Lee |239| . Concerted motion has also been implemented within lattice models 

1548 of exclusion processes in a more draconian manner. Gabel, Krapivsky, and Redner 

1549 have formulated a 'faciliated exclusion' process on a ring whereby a particle hops 

1550 forward to an empty site only if the site behind it is also occupied |241j . They find a 

1551 maximal current of 3 — 2-\/2 which is less than the maximal current of 1/4 arising in 

1552 the standard TASEP. A model that incorporates concerted motion might be described 

1553 by a facilitated exclusion model where motion in either direction occurs faster if the 

1554 particle is adjacent to another one. The hopping rate might also be increased if 

1555 an isolated particle moves to a site that results in it having an additional adjacent 

1556 particle. If these accelerated steps occur much faster than individual particle hops, 

1557 then the dynamics will resemble motion of pairs of particles. It would be interesting 

1558 to determine how this model of near-concerted hops increase the overall particle flux. 

1559 A second important point to emphasize is that different physical systems are 

1560 best modeled with varying degrees of asymmetry in the exclusion processes. Two 

1561 extreme limits are the totally asymmetric process, where particles only hop in one 

1562 direction, and symmetric exclusion, where particle hopping between any two adjacent 

1563 sites obeys detail-balance. In this case, detailed balance is violated only at the two 

1564 ends of the lattice. Net particle current arises only when there exists an asymmetry in 

1565 the injection and extraction rates at the two ends. This latter scenario corresponds to 

1566 boundary-driven exclusion processes where differences in the chemical potential of the 

1567 particles in the two reservoirs drive the flux. Osmotic and pressure-driven flows (when 

1568 local equilibrium thermodynamics holds and particle inertia is negligible) are examples 

1569 of processes best described using symmetric exclusion processes |242j . However, when 

1570 particles in the lattice are charged, and an external field is applied, the internal jumps 

1571 are asymmetric since there is a direct force acting on the particles, breaking detailed 
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1572 balance. A heuristic delineation between asymmetric and symmetric (boundary- 

1573 driven) exclusion processes can be motivated by considering the single-particle free 

1574 energy profiles. Figure [TU] depicts two hypothetical single-particle free energy profiles 

1575 experienced by an isolated particle under local thermodynamic equilibrium. 
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Figure 10. Representative single-particle, one-dimensional free energy 
landscapes over which excluding particles travel. The black curve free energy 
profile could represent a symmetric exclusion process, while the red landscape 
would describe an asymmetric exclusion processes. In the former, detailed balance 
is assumed to be violated at the entrance and exit sites and the particle flux is 
driven by differences between the chemical potentials of the two reservoirs. 

In biological systems, channels that transport water and neutral (e.g., sugars) or 
charge-screened molecules, can be described by symmetric exclusion, while motion of 
charged ions in single-file channels would be best characterized by partially asymmetric 
exclusions. Modeling uncharged particles allows equal internal hopping rates {p = q) 
with drive originating only at the boundaries The steady-state current is easily found 
and is a function of the asymmetry in the boundary injection and extraction rates 

" (i - l)(a + 5){I3 7) + p(a + /3 + 7 + ^) ■ ^^^^ 



The boundary injection and extraction rates are those defined in figure [TJc). As 
expected, this expression for J is similar to that arising from simple diffusive flux in the 
L oo limit. However, as discussed, nontrivial current fluctuations, which cannot be 

1587 accounted for by simple diffusion, have also been worked out |1611 11621 11641 12441 1245| . 

1588 When particles are 'charged', the level of asymmetry (the relative difference 
between an ion hopping forward and hopping backwards) is controlled by the 
magnitude of the externally applied 'electric field'. Besides the exact solutions for the 
steady-state current [125] , cumulants of the current in a weakly asymmetric exclusion 
process have also been derived |165| . the slowest dynamic relaxation mode computed 

1593 [1 71J , and the tracer diffusivity derived '163 . 

Specific biological realisations of driven transport include ion transport across ion 
channels [242) , while transport across nuclear pore complexes [2461 1247] , and osmosis 
[248 are typically boundary-driven, or ratcheted (see the section on molecular motors). 

In 1998, the X-ray crystal structure of the K+ ion channel was published Lj, 
allowing a more detailed mechanistic understanding of the ion conduction and 
selectivity across a wide class of related, and biologically important, ion channels |249) . 



+ This achievement was partially responsible for Peter Agre and Roderick McKinnonin being awarded 
the Nobel Prize in chemistry in 2003, "for discoveries concerning channels in cell membranes." 
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The crystal structure showed three approximately in-line vestibules that can each hold 
one K+ ion. Representations of the voltage-gated Kvl.2 channel are shown in figure 
111! Ion conduction can thus be modeled by a simple three-site partially asymmetric 
exclusion process provided the hopping rates at each occupation configuration are 
physically motivated [2481 12501 1251) . However, note that there is evidence that ion 
dynamics within some ion channels are 'concerted' because the fi'ee energy barriers 
between ion binding sites are small and ions entering an occupied channel can knock 
the terminal ion out of the channel. This 'knock-on' mechanism |252j has been 
motivated by molecular dynamics simulations |253l 1254] and studied theoretically via 
a reduced one-dimensional dynamical model |255) . 





(a) 



(b) 



Figure 11. (a) A ribbon figure of the Kvl.2 voltage-gated potassium ion channel 
embedded in a model lipid membrane. The pore structure clearly shows three 
dominant interior binding sites. This potassium channel image was made with 
VMD and is owned by the Theoretical and Computational Biophysics Group, 
NIH Resource for Macromolecular Modeling and Bioinformatics, at the Beckman 
Institute, University of Illinois at Urbana-Champaign. (b) A time series of a 
molecular dynamics simulation suggesting a concerted 'knock-on' mechanism. 
Image adapted from 254 with original labels removed, and used with their 
permission. 

1610 The symmetric exclusion process has also been used to consider osmotic flow 

1611 across membrane-spanning, single- file channels [248]. Here, the solvent flux is driven 

1612 by differences in the injection rates at the two pore ends connected to two separate 

1613 particle reservoirs. Simple osmosis across a strictly semipermeable membrane can be 

1614 described by a single species symmetric exclusion process where the injection rate of 

1615 the solvent from one of the reservoirs is reduced due to its smaller mole fraction arising 

1616 from the presence of solute in that reservoir |248j . The solvent current is approximately 

1617 linearly proportional to the difference in injection rates at the two ends of the lattice, 

1618 inversely proportional to the length of the lattice, but with a weak suppression due to 

1619 multiple pore occupancy |242| . 

1620 Another transport channel in cells are nuclear pore complexes (NPC), responsible 

1621 for the selective shuttling of large molecules and proteins (such as transcription factors, 

1622 histones, ribosomal subunits) through the double nuclear membrane. Exclusion 

1623 processes have been employed as theoretical frameworks for describing the efficiency 

1624 and selectivity NPC transport |247[|256) . 

1625 In all of the above applications, extensions exploiting multispecies exclusion 

1626 processes are often motivated. In biological systems ion transport is typically 'gated' 

1627 by cofactors that bind to the pore, either blocking ion transport, or inducing a 

1628 conformational change in the pore structure thereby affecting the entrance, exit, 

1629 and internal hopping rates |222| . The inclusion of additional species of particles in 

1630 the exclusion process has been used to describe 'transport factor' mediated nuclear 
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1631 pore transport |246[ 12471 1256) . For osmosis, solutes that are small enough to enter 

1632 membrane pores may also interfere with the transport of solvent particles through 

1633 channels. The solvent and solute species would have different entry, exit, and internal 

1634 hopping rates, describing their different interactions with the pore. In [248j . a simple 

1635 two-species, three-site, partially asymmetric exclusion model was used to show how 

1636 solutes that can enter a pore and suppress solvent flux. Moreover, it was found that a 

1637 small amount of slippage (passing of the solvent and solute particles over each other) 

1638 and a pore that was slightly permeable to solute can very effectively shunt osmotic flow. 

1639 When both solute and solvent can pass through the channel (with equal forward and 

1640 backward hopping rates in the interior), an interesting possibility arises whereby the 

1641 flux of one of the species can drive the other species from its low concentration reservoir 

1642 to the high concentration reservoir. Since the mechanism relies on entrainment of 

1643 particles that are driven up its chemical potential, slippage between the pumping 

1644 and convected particles destroys this entropic pumping mechanism. The efficiency of 

1645 using a pumping particle that travels from high chemical potential reservoir to low 

1646 chemical potential reservoir to pump the secondary particle was found using Monte- 

1647 Carlo simulations [257) . 




Figure 12. (a) Schematic of the Grotthuss proton transfer mechanism, (b) An 
MD simulation showing the persistence of a 'water wire' within a membrane- 
spanning Aquaporin-1 water channel protein [258) . The waters within the pore 
are shown magnified. 

1648 An even more interesting application of multispecies exclusion processes to 

1649 biological transport arises when the transported ion is a proton. In addition to 

1650 regulating ionic concentrations, pH regulation is an important aspect of cellular 

1651 function. It has been known for quite some time that the diffusivity of protons is 

1652 about an order of magnitude faster than that of small cations |2591 12601 . The physical 

1653 mechanism invoked to explain accelerated motion of protons is based on a proton 

1654 transfer or shuttling mechanism, analogous to electronic conduction. For a simple 

1655 ion to traverse a narrow ion channel, it must shed its hydration shell and push any 

1656 possible water molecules ahead of it within the pore. An extra proton, however, can 

1657 hop along an oxygen 'backbone' of a line of water molecules, transiently converting 

1658 each water molecule it visits into a hydronium ion HaO"*". The Grotthuss mechanism of 

1659 proton conduction has been conjectured to occur across many narrow pores, including 

1660 those in gramicidin- A channels |26H I262j , proton transfer enzymes such as carbonic 

1661 anhydrase j263j . voltage gated proton channels such as Hvl |264l . Aquaporin water 

1662 channels |265[ 12581 1266] , and carbon nanotubes (267, . All of these structures have in 
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1663 common the existence of a relatively stable water wire as shown on the right of figure 

1664 I12[ 



(a) 



ii. ^ H 
H H 



(b) 



l\ p- 

+ 1+1 O^+l + I - 



+ 

o 



+ 



-o 
o + 



o 



(c) 



o 


+ 






o 






o 






o 


o 




at 


t 




/ 




/ 










+ 






o 


+ 




o 


+ 




o 


+ 


o 








/ 


















tp 


+ 




o 


+ 


+ 




o 


+ 




o 


+ 


o 


o 




p+/ 


k 


-t 




\P- 






-t 




ty 


+ 





+ 


+ 








o 




o 




o 




St 










k 














o 





+ 


+ 






+ 


o 




o 


o 


+ 





Figure 13. (a) A cartoon of water- wire proton conduction. Hydrogen atoms and 
lone electron pairs of the water oxygen are shown, (b) A lattice model for the 
Grotthuss proton conduction mechanism. The symbols 0>+i s-nd ~ correspond 
to hydronium ions, water with lone pair pointing to the left, and water with 
lone electron pair pointing to the right, respectively. Proton movement from 
left to right leaves the donor site in a — configuration and can occur only if 
the receptor site is originally in the + configuration. The spontaneous left-right 
water flipping rates are k± and the proton forward and backward hopping rates 
(assuming a compatible configuration) are p±. (c) A time-sequence of a trajectory 
of configurations. Figures adapted from |268II269] . 

1665 The basic water wire mechanism can be mapped onto a three-species exclusion 

1666 model as shown in figure [T5i*l Unlike models where only one extra proton is allowed 

1667 in the channel at any given time |270| . the general exclusion model allows each site 

1668 to be in one of three states, protonated, water lone pair electrons pointing to the 

1669 left and to the right. In this particular application, no exact results are known, but 

1670 mean-field treatments have been used to extract qualitative phenomena (2681 1269) . 

1671 For example, in order to sustain a steady-state proton current, the lone-pair electrons 

1672 of a water molecule need to flip in order to relay successive protons. Not only was 

1673 proton conduction found to be mediated by water flipping, but nonlinear effects such 

1674 as saturation at high voltages and even negative differential resistances were exhibited 

1675 by the model [25511^ . 

* We caution the reader that the usage of 'M species' in the literature is not uniform or unique. 
Thus, 'three-species' here refer to three types of particles on the lattice, without holes. Thus, in 
terms of states at each site, it is the same as the 'two species' model describing solvent, solute, and 
holes at each site. 
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1676 The issue of concerted motions and precise definitions of the proton carrier species 

1677 also arises in more detailed considerations of proton transport. There has been 

1678 considerable effort devoted to identifying the precise molecular species that solvates 

1679 and relays the protons |271| . Moreover, there is some evidence that proton dynamics 

1680 in water wire conduction may be concerted [263]. Nonetheless, the three-species 

1681 exclusion model and its potential extensions are fruitful ways of understanding the 

1682 gross mechanisms in proton conduction. 

1683 4^.2. Simple models of molecular motors 

1684 Perhaps the simplest models of active biological transport are those of isolated 

1685 molecular motors that move along one-dimensional tracks such as actin, microtubules, 

1686 DNA, or RNA. Motors are enzymes such as dynein, kinesin, and myosin that hydrolyze 

1687 molecules such as ATP or GTP, turning the free energy released to directed motion 

1688 along their one-dimensional substrate (272[ 1273] . Motoring is necessary for sustaining 

1689 cell functions such as mediating cell swimming and motility and intracellular transport 

1690 of biomolecules, particularly at length scales where diffusion is not efficient, or 

1691 where spatial specificity is required. Due to the importance of molecular motors in 

1692 intracellular transport and cell motility, there is an enormous literature on the detailed 

1693 structure and chemo-mechanical transduction mechanisms of molecular motors. 

1694 From the theoretical physics point of view, molecular motors are useful examples 

1695 of non-equilibrium systems. Indeed, such motors typically operate far from 

1696 equilibrium, in a regime where the usual thermodynamic laws do not apply. A 

1697 molecular motor is kept far from equilibrium by a coupling with some external 

1698 'agent' (e.g., a chemical reaction, or a load-force): under certain conditions, work 

1699 can be extracted, although the motor operates in a medium with constant (body) 

1700 temperature. We emphasize that there is no contradiction with thermodynamics: the 

1701 system is far from equilibrium and the motor simply plays the role of a transducer 

1702 between the energy put in by the agent (e.g., chemical energy) and the mechanical work 

1703 extracted. Molecular motors have been described theoretically cither by continuous 

1704 ratchet models or by models based on master equations on a discrete space, which are 

1705 similar to exclusion-type systems, to which some of the exact analytical techniques 

1706 described above can be applied. Using an extension of the discrete two-state motor 

1707 model, Lau et al. [274] investigated theoretically the violations of Einstein and 

1708 Onsager relations and calculated the efficiency for a single processive motor in 

1709 [274 . Furthermore, it can be shown that the fluctuation relations (such as the 

1710 Gallavotti- Cohen theorem, the Jarzysnki-Crooks relations) play the role of a general 

1711 organizing principle. Indeed, cellular motors are systems of molecular size which 

1712 operate with a small number of molecules, and for these reasons undergo large thermal 

1713 fluctuations. The fluctuation relations impose general constraints on the function of 

1714 these nanomachines that go beyond classical thermodynamics. They provide a way 

1715 to better understand the non-equilibrium energetics of molecular motors and to map 

1716 out various operating regimes |274[ 12751 1276] . 

1717 Significant effort in the investigation into molecular motors has focussed on 

1718 identifying and understanding the molecular mechanics and the coupling of molecular 

1719 motion with a chemical reaction such as ATP hydrolysis. Form these studies, complex 

1720 descriptions of molecular motors have been developed, including a somewhat artificial 

1721 classification of motors employing 'power stroke' or 'Brownian ratchet' mechanisms. 

1722 This classification refers to how detailed balance is violated within the large motor 
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1723 molecule or enzyme. If certain internal degrees of freedom in a molecule are made 

1724 inaccessible at the right times, a net flux along these states can arise, ratcheting the 

1725 motion. If a particular transition is strongly coupled to specific chemical step such 

1726 as the hydrolysis of an ATP molecule bound in a pocket of the motor molecule, the 

1727 motion has been described as a power stroke motor. It has been shown that this 

1728 distinction is quantitative, rather than qualitative |277) . From a stochastic processes 

1729 point of view, when the dynamics of the microscopic internal motor degrees of freedom 

1730 are represented by, e.g., a Markov chain, power stroke and ratchet mechanisms can be 

1731 distinguished by where detailed balance is violated in the cycle. For example, if the 

1732 relative amount of violation occurs duriong transitions across states that are directly 

1733 coupled to motion against a load, the motor will be "ratchet-like" . 

1734 If the absicca in figure [TU] represents sequential internal molecular states, a 

1735 purely ratchet mechanism is naturally represented by a state evolving along the flat 

1736 energy profile. In this steady-state probability flux, or molecular motion, 

1737 arises from the absorbing and emitting boundary conditions that ratchet the overall 

1738 motion. In this stochastic picture, the power-stroke/Brownian ratchet distinction is 

1739 mathematically recapitulated by state-space boundary conditions and by the relative 

1740 amount of convection through state-space. How much a motor utilizes a power-stroke 

1741 mechanisms would be described by the Peclet number within the framework of single- 

1742 particle convection-diffusion or Fokker-Planck type equations. For a more detailed 

1743 delineation of regimes of mechanical-chemical coupling, see the review by Astumian 

1744 1278 . 

1745 In many cellular contexts, molecular motors are crowded and interact with each 

1746 other and exhibit collective behavior. Examples include connected myosin motors 

1747 in muscle, multiple motors and motor types on cellular filaments and microtubules 

1748 (often carrying large cargo such as vesicles), and motors that process DNA and RNA. 

1749 To model such systems, the details of how an isolated motor generates force and 

1750 moves may be best subsumed into a single parameter representing the mean time 

1751 between successive motor displacements. The internal dynamics, whether power- 

1752 stroke or Brownian ratchet, moves the motor one step against a load or resistance 

1753 at random times. These times are drawn from a distribution of hopping times that is 

1754 determined by the underlying, internal stochastic process. 

1755 One can simplify the modeling of molecular motors by assuming that the motor 

1756 stepping time is exponentially distributed with an inverse mean that defines the 

1757 hopping rate. Each motor hops along a one-dimensional track and can exclude other 

1758 motors. Typically, concerted motions are also neglected in this application. That is, 

1759 a motor is not allowed to push another one in front of it, moving both motors ahead 

1760 simultaneously. In the extreme limit, multiple motors are coupled with, e.g., elastic 

1761 elements and have been extensively modeled using simple Fokker-Planck equations 

1762 that incorporate mechanical and thermal forces [2801 EHD EHSl [2831 [SHU [2M1 [286] . 

1763 Weaker interactions that do not bound the distance between motors can be modeled 

1764 using the excluding particle picture by assigning different rules for the hopping of two 

1765 adjacent particles, analogous to the facilitated diffusion models of Gabel, Krapivsky, 

1766 and Redner [241] . Given the assumptions discussed above, exclusion processes 

1767 can be directly used to model collections of motors moving on one-dimensional 

1768 tracks. Additional effects particular to applications in biomolecular motors include 

1769 detachment and attachment kinetics and different types of motors that move in 

1770 opposing directions. 

1771 Molecular motors 'walk' along filaments but have a finite 'processivity' since 
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Figure 14. A schematic of two types of molecular motors moving along a cellular 
microtubule (adapted from 279 ). In this picture, kinesin moves to the right while 
dynein moves to the left. Each motor can be attached to large cargoes such as 
vesicles or other filaments. Note that the microtubule is constructed of twisted 
lanes of repeated molecular subunits. 



1772 they can spontaneously detach. Thus, they have a distribution of run lengths. 

1773 Conversely, motor molecules diffusing in the 'bulk' cytoplasmic space can attach to 

1774 interior locations of the lattice. Detachment violates particle conservation on the 

1775 lattice and has been studied using mean field models, hydrodynamic approximations 

1776 [IHZl [2HH1, and Monte-Carlo simulations HHl [211 USS- The detachment and 

1777 possible reattachment of particles on a lattice have been extensively studied in the 

1778 context of gas adsorption isotherms or 'Langmuir kinetics' (see |291| and references 

1779 within). While many models of Langmuir isotherms exist, previous studies considered 

1780 only passive, undriven particles. In the new work combining Langmuir kinetics with 

1781 driven exclusion processes, analytic progress can be made by considering the infinite 

1782 site, continuum hydrodynamic limit. If attachment and/or detachment occurs, the 

1783 hopping rates must be rescaled by the number of lattice sites such that the rates wa.d 

1784 (cf. figure [5]) at each site are inversely proportional to the number of sites L. If this is 

1785 not done, particles would only occupy a small region near the injection end of a long 

1786 lattice. To arrive at a nontrivial structure, the attachment and detachment rates must 

1787 be decreased so that the cumulative probability of desorption is a length-independent 

1788 constant. 

1789 As we mentioned in section |31 continuum hydrodynamic equations allow more 

1790 complicated models to be approximately treated. This has been the case for one- 

1791 dimensional exclusion processes with Langmuir kinetics. In the steady-state limit, 

1792 Parmeggiani et al. find [219j 

|^^ + (2p-l)^ + r!A(l-p(:r))-17Dp(x)=0, (86) 

1794 where r2A,D = ^k,dL = llia,d/£ represents appropriately rescaled detachment and 

1795 attachment rates that in this simple model are independent of L. A detailed asymptotic 

1796 analysis of equation (j86p was performed and a phase diagram as a function of four 

1797 parameters (the injection and extraction rates at the ends of the lattice, and the 

1798 adsorption and desorption rates) was derived [2181 1219) . They find a rich phase 

1799 diagram with coexisting low and high-density phases separated by boundaries induced 

1800 by the Langmuir kinetics. Langmuir kinetics have also been investigated in the 
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1801 presence of bottlenecks |292) . Klumpp and Lipowsky |293j have considered asymmetric 

1802 exclusion in tube- like structures such as those seen in axons or in filapodia [2941 1295) . 

1803 In this simplified geometry, the diffusion of motors within the tube can be explicitly 

1804 modeled [293 . 

1805 Since motors are used to carry cargoes across the different regions within a cell, 

1806 they travel on directed filaments and microtubules. Different motors are used to 

1807 carry cargo in one direction versus the other. For example, kinesins travel along 

1808 microtubules in the '+' direction, while dynein travels in the '-' direction, as shown 

1809 in figure 1141 This scenario can be modeled using a two-species exclusion process 

1810 with Langmuir kinetics overtaking. Moreover, microtubules are composed of multiple 

1811 twisted filamentous tracks. Not only can motors travel in opposite directions on the 

1812 same track, but filaments may be oppositely-directed within confined spaces such 

1813 as axons. Motors traveling in opposite directions on the same filament or on a 

1814 nearby parallel filament can be modeled using coupled chains of exclusion processes. 

1815 Many groups have explored the dynamical properties of two-species and two-lane 

1816 asymmetric exclusion processes (as shown in figure [71) |296[ 12971 12981 [2991 1300) . With 

1817 the proper biophysical identification of the parameters, results from these studies of 

1818 interacting lattices should provide illuminating descriptions of more complex scenarios 

1819 of moleculer motor- mediated transport. 

1820 4^.3. mRNA translation and protein production 

1821 One very special case of interacting motors moving along a one-dimensional track 

1822 arises in cell biology. In all cells, proteins are synthesized by translation of messenger 

1823 RNA (mRNA) as schematised in figure 1151 Complex ribosome enzymes (shown in 

1824 figure [TSl b')') unidirectionally scan the mRNA polymer, reading triplets of nucleotides 

1825 and adding the corresponding amino acid to a growing polypeptide chain. Typically, 

1826 many ribosomes are simultaneously scanning different parts of the mRNA. An electron 

1827 micrograph of such a polysome is shown in figure [TSl c). 

1828 The ribosomes are actually comprised of two subunits, each made up mostly 

1829 of RNA and a few proteins. Not only do these 'ribozymes' catalyze the successive 

1830 addition of the specific amino acid as they scan the mRNA, their unidirectional 

1831 movement from the 5' end to the 3' end, (where the 5' denotes the end of the mRNA 

1832 that has the fifth carbon of the sugar ring of the ribose as the terminus) constitutes a 

1833 highly driven process. Therefore, each ribosome is also a molecular motor that rarely 

1834 backtracks as it moves forward. The fuel providing the free energy necessary for 

1835 codon recognition and unidirectional movement is supplied in part by the hydrolysis 

1836 of GTP |301[ . Quantitatively, mRNA translation is different from typical molecular 

1837 motors in that ribosome processivity is very high, allowing one to reasonably neglect 

1838 detachments except at the termination end. Moreover, there are no known issues with 

1839 'concerted motions' or 'facilitated exclusion'. If one ribosome prematurely pushes the 

1840 one ahead of it, one would expect many polypeptides to be improperly synthesised. 

1841 Thus, at first glance, the TASEP seems to be the perfect model for this translation 

1842 process, with the particles being ribosomes and a site being a codon - a triplet of 

1843 nucleotides. On closer examination, it is clear that protein production is much more 

1844 complicated. Many biophysical features relevant to mRNA translation are missing 

1845 from the basic TASEP. The desire to have more 'realistic' models of protein production 

1846 has motivated the development of various extensions to the basic TASEP. In this 

1847 subsection, we will discuss only two efforts to generalize TASEP: allowing particles of 



NESM: A paradigm and applications 



50 




Figure 15. (a) A cartoon of the 'central dogma' in biology where mRNA 
is synthesized from cellular DNA and transported to the cytoplasm. The 
cytoplasmic mRNAs are then translated by ribosomes into polypeptides which 
may then finally be folded and processed into functioning proteins, (b) Crystal 
structure of both subunits of bacterial ribosome. (c) An electron micrograph of 
multiple ribosomes translating a single mRNA. 



1848 size £ > 1 and incorporating inhomogeneous, site-dependent hopping rates. In each 

1849 case, we wih describe the ceh biology which motivates the modifications. 

1850 Ahhough the fundamental step size of ribosomes is a codon, the large size of each 

1851 ribosome (20nm for prokaryotic ribosomes and 25-30nm for eukaryotic ribosomes) 

1852 means that they each cover approximately 10 codons along the mRNA chain and 

1853 exclude each other at this range. Therefore, TASEPs comprised of extended objects, 

1854 occluding £ > 1 lattice sites have therefore been developed, dating back to the late 

1855 1960's [131 mi- Although exact results for such a generalized TASEP on a ring are 

1856 available |302| . the problem with open boundaries remains unsolved. Instead, various 

1857 MFT's have been successful at capturing many important properties for modeling 

1858 mRNA translation. Let us devote the next paragraphs to this system. 

1859 First, consider a ring of L sites filled with TV particles, and note that every 

1860 configuration with particles of extent £ can be matched with a configuration with N 

1861 point particles {£ = 1) on a ring of L — N{£—1) sites. This mapping is easily understood 

1862 by regarding a configuration, C, as clusters of adjoining vacancies followed by clusters 

1863 of adjoining particles (regardless of their sizes). Therefore, the stationary distribution, 

1864 P* (C), is again flat and independent of £ for all p, q. Of course, the sum of the particle 

1865 density [p = N/L) and the hole density (phoio) now satisfy £p + phoic = 1, but with p 

1866 lying in a limited interval [0, !/£]■ With P* oc 1, finding the probability of particle-hole 

1867 pairs is straightforward, resulting in an exact expression [205) for the current density 

1868 relation, J{p) = /ophoio/(p + Phoic — ^/L), in the TASEP case. To lowest order in 1/L, 

1869 the formula 

l-[£-l)p 

1871 was known as early as 1968 [131 E], and leads to a maximal current of (1 -I- ^/£)^'^ 

1872 associated with the optimal density of p = {£ + V^)"^. Note that, though J{p) is no 

1873 longer symmetric about p, the particle-hole symmetry {p ^ Phoic) is preserved. This 

1874 invariance is expected on physical grounds, as the current arises only from exchanges of 

1875 particle-hole pairs. Though the stationary measure is trivial and G(r), the expectation 
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1876 value of two particles separated by r sites, can be written formally as a sum of products 

1877 of binomials. Of course, period-^ structures are to be expected in G (r). Despite the 

1878 conceptual simplicity of this problem, remarkably intricate patterns emerge |303) . 

1879 especially near complete filling, L = Ni |304j . Such structures are completely absent 

1880 from the £ = 1 case (whether in periodic or open TASEP), showing us that, even in 

1881 seemingly trivial situations, statistics of extended objects can produce surprises. 

1882 Turning to the open boundary TASEP, we must first specify how a particle 

1883 enters/exits the lattice. One possibility is 'complete entry, incremental exit' [305) . 

1884 where a particle may enter completely provided the first i sites are empty, while it may 

1885 exit moving one site at a time. In the NESS, exact expressions for the current like (I16p 

1886 can still be written. When both i and i+i are within [1, L], wc have J = (CTi(l — cri_|_^)). 

1887 From the 'incremental exit' rule, we have J = ((Ti_^) = . . . = (ct^-i) = /3 (ctl), leading 

1888 to a simple profile next to the exit. The major challenge comes from the 'complete 

1889 entry' condition: J = ^(111=1(1 ~ '^k))- Since this problem can no longer be solved 

1890 exactly, many conclusions can only be drawn from simulation studies or imposing mean 

1891 field approximations. For the latter, the nai've replacement of averages of products of 

1892 a by products of (a) (e.g., {uiUi+t) — > {ai){ai+e)) leads to extremely poor predictions. 

1893 Gibbs, et.al. [13l [M] took into account some of the effects of exclusion at a distance 

1894 and approximated J = ((Ti(l — Oi+i)) by 

J = (88) 

Pi+i + Pi 

1896 where = 1 — X^fetfi+i Pk is effective hole density in the i sites following i. For the 

1897 ring, the profile is flat and pt — p, so that ([88]) reduces to (I87p . Supplemented with 

1898 the appropriate boundary equations, ()88|) can be regarded as a recursion relation for 

1899 the profile. The successes of this approach include predicting a phase diagram that is 

1900 the same as the one in figure [5fc), except that the boundaries of the maximal current 

1901 phase are now sX a,f3 — ^1 + • Simulations largely confirm such predictions 

1902 [2051 1305] . suggesting that the correlations neglected by the scheme of Gibbs, et al. 

1903 |13[ 114] are indeed small. On the other hand, for more sensitive quantities like the 

1904 profile, this MET is less successful, especially for the high density phase (/? ^ a ^ 1). 

1905 To produce better predictions, Shaw, et al. |103j introduced a more sophisticated 

1906 MET, taking into account some pair correlations. So far, no higher level of MET have 

1907 been attempted. 

1908 The second modification to the basic TASEP we consider here is site-dependent 

1909 hopping rates. Since the translation of mRNA into polypeptides depends on the 

1910 sequence of nucleotides, the hopping rates of the ribosome TASEP particles can vary 

1911 dramatically as a function of its position on the lattice. The local hopping rates 

1912 depend on the effective abundance of the specific amino-acid-charged tRNA that 

1913 participates in each elongation step at each site. One of the first treatments of 

1914 TASEPs on a nonuniform lattice considered a single defect, or slow hopping sited, 

1915 in the middle of an asymptotically long lattice. Kolomeisky derived the expected 

1916 steady-state particle flux across a lattice with a single slow hopping site by ignoring 

1917 particle-particle correlations across the defect [306j . He self-consistently matched the 

1918 exit rate from the flrst half of the chain with the injection rate of the second half of 

1919 the chain. This approach indicated that each long uniform region between slow sites 

tl If an isolated site had faster hopping than all its neighbors, the average current is hardly affected, 
though there are noticable changes to the profile. Thus, we focus on slow sites. 
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1920 can approximately be treated as separate, but self-consistently connected, uniform 

1921 exclusion processes. Later work by Chou and Lakatos developed a more refined method 

1922 of connecting the current between sections separated by interior defects. Their method 

1923 generalizes mean-field theory by explicitly enumerating the configurations of the sites 

1924 straddling the inhomogeneity, and self-consistently matching the current across this 

1925 segment with the asymptotically exact steady-state currents across the rest of the 

1926 homogeneous lattice. Steady-state particle currents across localized defects can be 

1927 accurately computed [307]. Using this approach, the synergistic effects of two nearby 

1928 slow hopping sites were analyzed. It was shown that two defects decreased the current 

1929 more when they are placed close to each other. The decrease in current approaches 

1930 that of a single defect as the distance between two defects diverges since the dynamical 

1931 'interaction' from overlap of density boundary layers vanishes. 




Figure 16. mRNA translation with fixed slow codons, or bottlenecks across 
which ribosomc motion is slower (p* < p) than across other codons. This local 
slowing down can arise from a limited supply of the appropriate amino acid- 
charged tRNA corresponding to the slow site. 

1932 Another method for approximately treating inhomogeneous hopping rates was 

1933 developed by Shaw et al. [103j where the original mean-field equations of Gibbs 

1934 et al. [131 [2] were generalized to include site-dependent elongation rates. Dong, 

1935 Schmittman, and Zia [308[ 13091 systematically analyzed the effects of defects on the 

1936 steady-state current of extended particles that occupy multiple (w > 1) sites. They 

1937 used a combination of self-consistent mean field theory and Monte-Carlo simulations 

1938 to show the importance of the location of the defect, especially if the defect is placed 

1939 close entry or exit sites where they may possibly 'interact' with the end-induced density 

1940 boundary layers. 

1941 As mentioned in section [5751 hvdrodvnamic MFT developed for analyzing systems 

1942 with slowly varying hopping rates. Continuum equations, such as (|84p . that 

1943 incorporate spatially varying hopping rate functions p{x) and q{x) have been derived 

1944 for the mean field ribosome density p{x) [2051 12061 12071 ITM] . For single-site particles, 

1945 these hydrodynamic equations can be integrated once to arrive at singular nonlinear 

1946 first order differential equation where the integration constant is the steady-state 

1947 particle current. Using singular perturbation theory, analytic solutions for density 

1948 profiles have been found for very specific hopping rate functions p{x) and q{x) 

1949 [124 1 1204 ) . 
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1950 The stochastic process underlying the models of niRNA translation described 

1951 above all assume an exponentially-distributed waiting time between elongation events. 

1952 Dwell time distributions would be difficult to incorporate into simple ODEs for 

1953 the mean-field particle density or continuum hydrodynamic equations. Although 

1954 different dwell time distributions would not qualitatively affect steady-state ribosome 

1955 current and protein production rate, to model them explicitly requires incorporation 

1956 of intermediate biochemical steps involved in initiation, elongation, and termination 

1957 of the individual ribosomes. These steps include binding of an amino-acid charged 

1958 transfer RNA (tRNA) to one of the binding sites, hydrolysis, etc. [3101 . Models 

1959 that include more complex elongation steps have also been developed by a number 

1960 of researchers [ml |3ll |313l [2261 [SH IMS |315]. The main qualitative result of 

1961 these studies is that the standard current phase diagram is shifted due to a varying 

1962 effective elongation (internal hopping) rate. Since the internal hopping rate depends 

1963 on the details of their models, it cannot be nondimensionalized and the standard phase 

1964 diagram is roughly reproduced with aoff/Pcff and Pcg/PcG as the tuning variable. Here 

1965 acff, /3off , and PcS are are effective rates that might be estimated as the inverse of the 

1966 mean of the associated dwell time distributions. 




(a) (b) 



Figure 17. (a) AFM image of circularized mRNA. (b) Schematic of a model 
where ends of mRNA are sticky (the poly-A tails are known to bind to initiation 
factors [3161 13171 1318| ). increasing the probability of loop formation. 

1967 Finally, the mRNA translation process, much like molecular motor-facilitated 

1968 intracellular transport, occurs in complex, spatially heterogeneous environments and 

1969 involves many molecular players that initiate and terminate the process. For example, 

1970 certain initiation factors that prime the initiation site for ribosomal entry have been 

1971 shown to bind to the polyadenylated tails of mRNA, thereby forming circularized 

1972 RNA. One hypothesis is that mRNA circularization facilitates the recycling of 

1973 ribosomes. In circularized mRNA, ribosomes that detach from the termination site 

1974 after completing translation have a shorter distance to diffuse before rebinding to 

1975 the initiation site. A model for this effect was developed in |319| . where an effective 

1976 injection rate was self-consistently found by balancing the steady-state ribosome 

1977 concentration at the initiation site with the diffusive ribosome fiux emanating from 

1978 the nearly termination site. In this model, the diffusive feedback tends to increase 

1979 the protein production rate, especially when overall ribosome concentrations are low. 

1980 One can also imagine a strong feedback if factors regulating translation initiation were 
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1981 
1982 
1983 
1984 
1985 



1990 
1991 
1992 
1993 
1994 



themselves products of translation. Newly produced cofactors can readily maintain 
the initiation rate a at a high value. 

An even more important aspect of mRNA translation is that ribosomes and 
initiation factors |320) . and mRNAs |321| can be actively localized to endoplasmic 
reticular (ER) membranes and compartments, depending on what types of proteins 
they code for, and where these proteins are needed. In confined cellular spaces, the 
supply of ribosomes and initiation factors may be limited. Moreover, there are many 
different mRNA copies that compete with each other for ribosomes. This global 
competition has been modeled by Cook et al. |322[ 1323] who defined an effective 
initiation rate aes{N) which is a monotonically increasing function of the free ribosome 
concentration. They considered mRNAs of different length (but of identical initiation, 
elongation, and termination rates) and found that steady-state protein production 
for different length niRNA's were comparable, but that their ribosome loading levels 
exhibit varying levels of sensitivity on the total ribosome mass. 



1995 4.4. Free boundary problems and filament length control 



1996 
1997 
1998 
1999 
2000 
2001 



Our final example of a class of biological application of exclusion processes involves 
changes in the length of the underlying one-dimensional substrate. A dynamically 
varying lattice length arises in at least three different cellular contexts: growth of 
filaments such as hyphae and cellular microtubules, replication forks, and mRNA 
secondary structure, in each of these examples, there is a 'moving boundary' whose 
dynamics depends on transport within the domain bounded by the boundary. 



(a) 
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Figure 18. (a) Schematic of multiple TASEPs a dynamically growing tips. 
The particles are motors carrying building blocks of the underlying filaments 
[228| . (b) A model of microtubule length control where particles reaching the tip 
depolymerises the last subunit 12311 . Figures are adapted from those in |228l and 
[2MI. 



An analogy can be made with the classic free boundary problem arising in 
continuum physics. In the 'Stefan problem' a diffusive quantity mediates the growth 
of an interface bounding the dynamics. This diffusing quantity may be heat which 
melts away water-ice interface [324], or a particulate systems that deposits and 
makes the interface grow [233j . In systems relevant to cell biology, the dynamics 
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2007 (i) occur in confined, often one-dimensional geometries, and (ii) occur on a fluctuating 

2008 mesoscopic or molecular scale. In confined geometries, particle-particle exclusion 

2009 become important. By contrast, this exclusion is absent when the transported field 

2010 is say, heat. On small scales, statistical fluctuations of the interface, and its coupling 

2011 with the fluctuating particle field can also be important. 

2012 A stochastic one-dimensional moving boundary problem can also be described 

2013 within an exclusion process where the number of sites is allowed to fluctuate. Such 

2014 models have been described in [2321 EMI ESU [2321 [230], and [231]. In [232], a 

2015 fluctuating wall is coupled to an asymmetric exclusion process as shown in figure 

2016 [9] The wall particle acts like a piston and defines one boundary of the lattice and 

2017 has its own intrinsic forward and backward hopping rates. Particles impinge on the 

2018 wall and have a certain detachment rate. Using a moving frame version of asymptotic 

2019 matching [307] the expected position and fluctuations of the wall were accurately 

2020 computed. The wall was found to either extend the lattice indefinitely, compress the 

2021 particles and fall off the lattice at the injection site, or find an equilibrium lattice 

2022 length. 

2023 In terms of specific biological applications, Evans and Sugden |228j and Sugden 

2024 et al. |229j consider a free boundary TASEP as a model for hyphae growth in 

2025 fungi (cf. figure ITSf a)). The hypothesis is that kinesin motors carry cargo and 

2026 move unidirectionally toward the tip. The delivered cargo can then extend the tip 

2027 by incrementing the number of lattice sites of the TASEP. In this case, there is no 

2028 confining wall, and the velocity tip extension is proportional to the number of arriving 

2029 TASEP particles. Even though the tip is always growing with a fixed velocity, the 

2030 particle density profiles can acquire different structure, including a jammed, high 

2031 density phase. 

2032 Another application of dynamic-boundary TASEPs involving molecular motors 

2033 was described by Hough et al. |231j . In their model of microtubule length control, 

2034 kinesin-8 motors move along a microtubule and depolymerize the tip when it is reached. 

2035 The tip also has an intrinsic degradation and assembly rate, the latter depending on 

2036 the concentration of subunits in the bulk. Therefore, this model is similar to a TASEP 

2037 with a confining wall |232| representing the end of the lattice. However, in this model, 

2038 the kinesis motors also undergo Langmuir kinetics by attaching to, and detaching 

2039 from the microtubule lattice. Hough et al. find regimes where microtubule length- 

2040 dependent depolymerization rates arise, as well as how the behavior depends on bulk 

2041 motor concentrations [231] . 

2042 A number of potentially new models remain relatively unexplored. For example, 

2043 molecular motors often encounter an obstacle. As the helicase motor separates DNA 

2044 strands for transcription, it must break base pair bonds and push the replication fork 

2045 forward }325[ I326| I327| I328j . The fork may be represented by a confining wall that 

2046 tends to reseal the single strands of DNA. Similarly, mRNA translation by ribosome 

2047 'motors' typically occur in the presence of hairpins that must be separated before the 

2048 ribosomes can progress. In both of these applications the motor-wall interaction may 

2049 be considered passive (ratchet-like) or active (forced separation) [328 . These limits 

2050 are again related to the issue of concerted motion or faciliated exclusion mentioned 

2051 earlier. The motor can actively advance the replication fork, whereby the forward 

2052 motions of the leading particle and the wall occur simultaneously in a concerted 

2053 fashion. Alternatively, the opening of the fork can be thermally activated, temporarily 

2054 allowing the leading motor to ratchet the wall. Sometimes, due to the particular 

2055 geometry, a motor may not need to move against a barrier, but may pass through it. 
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2056 An interesting realisation arises in mRNA translation where ribosome may jump over 

2057 hairpins, effectively making hops to an empty site far away on the lattice. This kind 

2058 of 'short-cut' model has been studied using mean-field theory and simulations by Kim 

2059 et al. |329| and have been found to exhibit coexistence of empty and jammed phases. 

2060 The consequences of these rich mechanisms on cell biology remain to be explored. 

2061 5. Summary and outlook 

2062 Non-equilibrium statistical mechanics is the study of stochastic processes for a system 

2063 involving many interacting degrees of freedom. As such, it concerns essentially the 

2064 entire spectrum of natural phenomena, ranging over all length scales and relevant for 

2065 all areas in science and engineering. While a small encyclopedia would be needed to 

2066 adequately cover NESM, reviews dealing with more specialised topics are abundant 

2067 in the literature. In this article, we attempt a compromise by providing (a) a bird's 

2068 eye view of NESM, and (b) a 'bridge' between the two different approaches to the 

2069 subject. For the former, we review some of the fundamental issues associated with 

2070 NESM and why the techniques (and assumptions) of equilibrium statistical mechanics 

2071 cannot be easily generalized. Using the language of stochastic processes, we base our 

2072 discussions on master equations and focus on systems which violate detailed balance 

2073 (or time reversal). One serious consequence is that, given a set of such rates, even 

2074 the stationary distribution, P* - a NESS analog of the Boltzmann factor e~^^ - 

2075 is generally unknown. There is an additional challenge: the presence of persistent 

2076 probability currents (forming closed loops, of course), K*, much like the steady electric 

2077 current loops in mangetostatics. By contrast, K* = if the transition rates obey 

2078 detailed balance, in analog to electrostatics. We emphasized that these currents play 

2079 important roles in a NESS, as they lead to persistent fluxes of physically observable 

2080 quantitities (e.g., mass, energy, etc.) as well as constant entropy production in the 

2081 surrounding medium. 

2082 Naturally, our 'bird's eye view' is both limited and colored. For example, we 

2083 offer only a brief glance on fluctuation theorems and symmetry relations, a topic of 

2084 NESM which has witnessed considerable progress over the past two decades. Even 

2085 more generally, our focus on master equations implicitly assumes that the dynamics 

2086 can be decomposed into a series of Markov processes with exponentially distributed 

2087 waiting times in each conflguration. Age-dependent processes described by more 

2088 complex waiting time distributions of each configuration may often be better studied 

2089 using theories of branching processes. Moreover, master equation approaches are 

2090 not very suited for studying the statistics of individual trajectories or realisations. 

2091 Understanding properties of individual trajectories in configuration space is important, 

2092 for example, in problems of inference and recontruction of the stochastic rules from 

2093 data. These attributes are better analysed using techniques of, e.g., stochastic calculus 

2094 or path integrals. Therefore, a comprehensive overview of the foundations of NESM 

2095 and stochastic processes would include concepts and formalism from other disciplines 

2096 such as probability theory, statistics, and control theory - a task far beyond our scope. 

2097 To meet the challenges posed by NESM, a wide range of approaches have been 

2098 developed. On the one hand, we have model systems, sufficiently simple that rigorous 

2099 mathematical analysis is feasible. At the other extreme are models which account for 

2100 many ingredients deemed essential for characterizing the large variety of phenomena 

2101 found in nature. Understandably, those working on either end of this spectrum are not 

2102 typically familiar with the progress at the other end. A historic example appeared in 
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2103 the late 1960's when Spitzer [5] and Gibbs, et.al. [131 [13] investigated independently 

2104 the same stochastic process, but were unaware of each others' studies. In this article, 

2105 we attempt to bridge this divide by providing brief reviews of a limited set of topics 

2106 from both ends. 

2107 For the exclusion process, we presented many exact results, along with two 

2108 complementary techniques to obtain them. The stationary distribution for a system 

2109 with periodic boundary conditions (i.e., ASEP on a ring) was known to be trivial [9]. 

2110 On the other hand, a rich variety of dynamic properties, even for SEP, have been 

2111 discovered over the last four decades. With open boundaries, particles may enter 

2112 or exit the lattice at both ends, so that the total occupation becomes a dynamic 

2113 variable. As a result, this system poses much more of a challenge. Little was known 

2114 until the 1990's, when the exact P* was found |118[ I120[ 1121) through an ingenious 

2115 recognition of a matrix representation. For simplicity, we presented the details of 

2116 this method only for the extreme case of TASEP, in which particles enter/exit with 

2117 rate a/ (3 (section [3. 4p . Once P* is known, it is straightforward to compute the exact 

2118 particle current, J, and the mean density profile pi. A non-trivial phase diagram 

2119 in the a-/3 plane emerged, with three distinct phases separated by both continuous 

2120 and discontinuous transitions (section 13. 5|) . In addition, algebraic singularities are 

2121 always present in one of the phases. Such phenomena are novel indeed, since the 

2122 conventional wisdom from EQSM predicts no phase transitions (i.e., long-range order), 

2123 let alone generic algebraic singularities, for systems in one dimension, with short- 

2124 ranged interactions! Beyond NESS, another powerful method - the Bethe Ansatz 

2125 ~ have been applied successfully to obtain time-dependent properties f section l3.7p . 

2126 Analyzing the spectrum and the eigenvectors (Bethe wave functions) of the full Markov 

2127 matrix, other physically interesting quantities in a NESS can be also computed. 

2128 Examples include the fluctuations of the current (around the mean J), encoded in 

2129 the large deviation function, G{j). The knowledge of dynamic properties also allows 

2130 us to explore, in principle, ASEPs on infinite lattices, for which the initial configuration 

2131 must be specified and the system properties after long times are of interest. In practice, 

2132 however, these problems are attacked by another technique, fascinatingly related to the 

2133 theory of random matrices, representations and combinatorics (section 13.101) . These 

2134 powerful methods have yielded mathematically elegant and physically comprehensible 

2135 results, rendering exclusion models extremely appealing for modeling real systems. 

2136 Apart from exact solutions, we also presented an important approximation scheme 

2137 ~ the mean-field theory. Based on physically reasonable and intuitive arguments, 

2138 we consider coarse-grained densities as continuum fields, obeying hydrodynamic 

2139 equations, e.g., the density p (a;, t) in ASEP. Remarkably, such a MET predicted some 

2140 aspects of TASEP exactly (e.g., the current density relation J = p(l — p)). Though 

2141 neither rigorous nor a systematic expansion (such as perturbation theory with a small 

2142 parameter), MET has provided valuable insights into behaviour which we cannot 

2143 compute exactly. This is especially true when appropriate noise terms are added, 

2144 so that we can access not only the (possibly inhomogeneous) stationary profiles but 

2145 also the fluctuations and correlations near them. A good example is the average power 

2146 spectrum associated with N{t), the total mass in an open TASEP. Intimately related 

2147 to the time-dependent correlation of the currents at the two ends, this quantity has 

2148 eluded efforts to find an exact expression, despite our extensive knowledge about the 

2149 full Markov matrix. Instead, starting with a MET and expanding around a fiat profile 

2150 in the high/low density phase, many interesting properties of this power spectrum 

2151 can be reasonably well understood |330| 1331) . Undoubtedly, there are many other 
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2152 physically interesting quantities (in ASEP and other 'exactly solvable systems') for 

2153 which theoretical predictions can be obtained only from mean field approximations. 

2154 Of course, as we consider modeling transport phenomena in biological systems, more 

2155 realistic features must be included. Then, MFT typically offers the only route towards 

2156 some quantitative understanding. 

2157 We next turned to specific generalizations of the ASEP which are motivated by 

2158 problems in biological transport. While the exactly solvable ASEP consists of a single 

2159 lattice of fixed length with at most one particle occupying each site and hopping 

2160 along from site to site with the same rate, the variety of transport phenomena in cell 

2161 biology calls for different ways to relax these constraints. Thus, in the first model for 

2162 translation in protein synthesis, Gibbs, et.al. already introduced two features 114] 

2163 absent from the standard ASEP. Representing ribosomes/codons by particles/sites, 

2164 and being aware that ribosomes are large molecules covering 0(10) codons, they 

2165 considered ASEP with extended objects. They were also mindful that the sequence of 

2166 codons in a typical mRNA consists of different codons, leading to a ribosome elongating 

2167 (moving along the mRNA) with possibly different rates. The result is an ASEP 

2168 with inhomogeneous hopping rates. Beyond these two aspects, we know that there 

2169 are typically thousands of copies of many different mRNA's (synthesizing different 

2170 proteins) within a single cell. Now, they compete for the same pool of ribosomes. 

2171 To account for such competition, we should study multiple TASEPs, with different 

2172 lengths, 'interacting' with each other through a common pool of particles. Within 

2173 the cell, the ribosomes presumably diffuse around, leading to possibly more complex 

2174 pathways of this 'recycling of ribosomes.' 

2175 Of course, ribosome motion along mRNA is only one specific example of molecular 

2176 motors. Therefore, the TASEP is also suited for modeling cargo transport by molecular 

2177 motors 'walking' along cellular microtubules. In this case, many essential biological 

2178 features are also absent from the standard TASEP model. Motors are known to 

2179 detach from the microtubule and reattach at other points, so that Langmuir kinetics 

2180 [291] should be introduced. The cargoes they carry are typically much larger than 

2181 their step size, leading us again to long-ranged particle-particle interactions. There 

2182 are many lanes on a microtubule, so that we should include multiple TASEPs in 

2183 parallel, with particles transfering from one lattice to another, much like vehicular 

2184 traffic on a multilane highway. Molecular motors come in many varieties (dynein, 

2185 kinesin, etc.) which move in opposite directions and at different speeds. Consequently, 

2186 a proper model would consist of several particle species hopping in different directions. 

2187 The variety of speeds is the result of complex, multi-step chemical reactions, so that 

2188 the dwell times are not necessarily distributed according to simple exponentials. To 

2189 account for such details, particles with internal states can be used necessary. This 

2190 level of complexity is also sufficient for modeling another important biological process: 

2191 transport through channels on membranes (pores). Various ions, atoms and molecules 

2192 are driven in both directions, often 'squeezing by' each other. Finally, microtubles grow 

2193 and shrink, a process modeled by a dynamic L, the length of our lattice. Typically, 

2194 the associated rates are governed by the densities of particles at the tip, leading us to 

2195 an entirely new dimension in mathematical complexity. 

2196 Finally, let us provide an outlook of NESM beyond the topics presented in the 

2197 sections here. In the realm of exactly solvable models, ASEP is just one of many. 

2198 Not surprisingly, all but a few are one-dimensional systems. For example, we noted in 

2199 section [3791 the zero-range process. Also introduced by Spitzer |9], it is a closely related 

2200 model for mass transport. Multiple occupancy is allowed at each site, while some of the 
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2201 particles hop to the next site. Thus, this process is well suited to describe passengers 

2202 at bus stops, with some of them being transported to the next stop. Especially 

2203 interesting is the presence of a 'condensation transition,' where a macroscopic number 

2204 of particles occupy a single site as the overall particle density on a ring is increased 

2205 beyond a critical value. Much progress has emerged, especially in the last two decades 

2206 (see |332j for a review). Another notable example of transport is the ABC model 

2207 IHH [87], mentioned at the end of section [2] In general, it also displays transitions 

2208 of a non-equilibrium nature, admitting long-range order despite evolving with only 

2209 short-ranged dynamics. 

2210 Apart from transport models, exact results are also known for systems with no 

2211 conservation law. A good example is the kinetic Ising chain [53], coupled to two 

2212 thermal reservoirs at different temperatures (T > T'). As a result, it settles down 

2213 to a NESS, with generally unknown P*. Depending on the details of the coupling, 

2214 exact results are nevertheless available. In particular, if every other spin is updated 

2215 by a heat-bath dynamics associated with T and T', then all multi-spin correlations 

2216 are known exactly |3331 1334] . Remarkably, even the full time dependence of these 

2217 correlations can be obtained exactly, so that the full P(C, t) can be displayed explicitly 

2218 as well [335] IB- As a result, it is possible to compute the energy flow through the 

2219 system (from the hotter bath to the cooler one) in the NESS, as well as the entropy 

2220 production associated with the two baths |333| . Similar exact results are available in 

2221 a more common form of imposing two baths, namely, joining two semi-finite chains 

2222 (coupled to T and T') at a single point. Since the energy flows from the hotter bath, 

2223 across this junction, to the cooler side, we may ask for say, (a) the power injected to 

2224 the latter and (b) the profile of how this injected energy is lost to the colder bath. Only 

2225 the average of the latter is known [337j , but the large deviation function of the total 

2226 injected power can be computed exactly |338[ 1339] . These are just some examples of 

2227 other exactly solvable systems which evolve with detailed balance violating dynamics. 

2228 In the realm of potential applications, exclusion processes extend well beyond the 

2229 examples in biological transport presented here. In the general area of 'soft condensed 

2230 matter,' the exclusion mechanism arises in many other systems, such as motion of 

2231 confined colloids |3401I341] . Further afield, the process of surface growth in a particular 

2232 two-dimensional system can be mapped to an ASEP [ 342i 1182) . On larger scales, 

2233 exclusion processes lends themselves naturally to modeling traffic flow [3431 13441 1345] 

2234 and service queues [346] . For each of these applications, though ASEP and its variants 

2235 may not be sufficiently 'realistic', they nonetheless provide a succinct physical picture, 

2236 some insights from mean-field analysis and a precise language on which sophisticated 

2237 mathematical techniques can be applied. Along with the overall improvement of 

2238 various technologies (in e.g., nanoscience, renewable energy), we expect that there 

2239 will be many opportunities for exclusion processes to play a role, both in modeling 

2240 newly discovered phenomena and in shaping directions of further research. 

2241 Broadening our outlook from exactly solvable models and potential applications 

2242 to NESM in general, the vistas are expansive. It is beyond our scope to provide 

2243 an exhaustive list of such systems, which would include problems in aging and 

2244 branching processes, directed percolation, dynamic networks, earthquake prediction, 

2245 epidemics spreading, granular materials, financial markets, persistence phenomena, 

2246 population dynamics, reaction diffusion, self organized critically, etc. On the purely 

ft Indeed, due to the simplicity of heat-bath dynamics, some exact results are known even if the rates 
are time dependent |336j . 
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2247 theoretical front, many fundamental issues await further exploration. For example, 

2248 the implications of probability currents, beyond the computation of physical fluxes, 

2249 may be far reaching. If we pursue the analog with electromagnetism, we could ask if 

2250 these currents can be linked to a form of 'magnetic fields,' if there is an underlying 

2251 gauge theory, and if these concepts arc constructive. Perhaps these ideas will lead 

2252 us to a meaningful framework for all stationary states, characterized by the pair 

2253 of distributions {P*,K*}, which encompasses the very successful Boltzmann-Gibbs 

2254 picture for the equilibrium states. In particular, in attempting to describe systems 

2255 which affect, and are affected by, their environment (through e.g., entropy production) 

2256 NESS represents a significant increase of complexity from EQSS. Of course, the hope 

2257 is that an overarching theory for NESM, from the full dynamics to predicting NESS 

2258 from a given set of rates, will emerge in the near future. Such a theory should help 

2259 us reach the ultimate goal for say, biology - which would be the ability to predict the 

2260 form and behaviour of a living organism, based only on its intrinsic DNA sequence and 

2261 the external environment it finds itself. For the latter, we have in mind both sources 

2262 of input (e.g., light, air, food, stimulations) and output (e.g., waste disposal, work, 

2263 reproduction). In the absence of such interactions with the environment, an isolated 

2264 DNA will evolve to an equilibrium state - probably just an inert macromolecule. 

2265 To fully understand the physics of life, we believe, a firm grasp of non-equilibrium 

2266 statistical mechanics is absolutely vital. 
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